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Preface 


This  report  is  the  final,  summary  report  for  Work  Unit  CU-1062,  “Develop¬ 
ment  of  Simulators  for  In  Situ  Remediation  Evaluation,  Design,  and  Operation,” 
which  was  under  the  Cleanup  Pillar  of  the  Strategic  Environmental  Research  and 
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Unit  CU-1062  was  conducted  by  U.S.  Army  Engineer  Research  and  Development 
Center  (ERDC),  imder  the  purview  of  the  Enviromnental  Laboratory  (EL), 
Vicksburg.  MS.  ERDC  involvement  in  SERDP  was  coordinated  by  Dr.  John 
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Ms.  Catherine  Vogel  was  the  Cleanup  Program  Manager  for  SERDP. 

The  technical  leader  and  point  of  contact  for  Work  Unit  CU-1062  was 
Dr.  Mark  S.  Dortch,  Chief,  Water  Quality  and  Contaminant  Modeling  Branch 
(WQCMB),  Environmental  Processes  and  Engineering  Division  (EPED),  EL. 
Other  principal  investigators  participating  in  the  project  included  Mr.  Christian  J. 
McGrath,  WQCMB;  Dr.  Carlos  E.  Ruiz,  WQCMB;  Dr.  Stacy  E.  Howington, 
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University,  Blacksburg;  University  of  Alabama,  Tuscaloosa;  and  Tennessee 
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1  Introduction 


Background 

The  Departments  of  Defense  (DoD)  and  Energy  (DOE)  are  involved  in  the 
cleanup  of  Aousands  of  sites  at  an  estimated  cost  of  over  hundreds  of  billions  of 
dollars.  Significant  resources  will  be  consximed  in  the  evaluation  of  alternative  or 
combined  remediation  technologies,  particularly  where  in  situ  or  innovative 
technologies  are  considered.  Yet  there  is  no  expedient  means  to  quantitatively 
evaluate  alternatives.  Such  methods  could  provide  more  accurate  and  technically 
defensible  evaluations.  Expensive  remediation  decisions  based  on  inadequate 
information  will  become  even  more  expensive  if  they  result  in  inefficient  design 
or  failure.  These  failures  may  result  in  increased  human  and  environmental  risk, 
thereby  reducing  public  confidence  and  negatively  impacting  DoD  and  DOE 
readiness.  Methods  to  streamline  the  evaluation  process  would  also  save  time  and 
money.  The  ultimate  goal  in  remediation  modeling  is  to  minimize  remediation 
costs  and  environmental  and  human  health  risks  while  maximizing  cleanup. 
Toward  this  end,  the  general  goals  of  this  project  were  to  (a)  develop  reliable 
simulators  for  promising  technologies  of  interest  to  DoD,  DOE,  and  the  regulatory 
community  including  the  U.S.  Environmental  Protection  Agency  (EPA);  and 
(b)  provide  efficient  access  to  multiple  remediation  simulators  through  a  common 
user  environment  amenable  to  multidisciplinary  cleanup  teams.  A  common 
graphical  user  environment  has  been  developed  for  these  simulators:  the  DoD 
Grovmdwater  Modeling  System  (GMS). 

The  GMS  provides  conceptualization,  characterization,  model  parameteri¬ 
zation  and  setup,  visualization,  and  animation  capabilities.  GMS  extensions,  either 
ongoing  or  planned,  will  provide  capabilities  for  conducting  parameter  estimation, 
rmcertainty  analysis,  optimization,  and  cost  analyses.  The  U.S.  Army  Engineer 
Research  and  Development  Center  (ERDC)  leads  the  GMS  development,  which 
has  been  a  Tri-Service,  multiagency  activity  involving  multiple  DoD,  DOE,  and 
EPA  researchers  and  20  university  partners.  The  GMS  is  employed  by  himdreds 
of  users  within  the  Tri-Services,  DOE,  and  EPA.  Commercial  versions  of  the 
GMS  are  used  by  many  groups  within  the  cleanup  industry. 
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Objectives 

The  overall  objective  of  this  project  was  to  provide  DoD,  DOE,  EPA,  and 
other  GMS  users  with  the  ability  to  simulate  in  situ  remediation  and  natural 
attenuation  of  subsurface  contamination.  This  capability  allows  the  evaluation  of 
the  effectiveness  of  various  cleanup  strategies  and  improvement  of  system  design 
prior  to  implementation.  Simulators  also  provide  a  means  to  evaluate  cleanup 
system  operational  strategies  and  improvements.  Additionally,  products  from  this 
project  provide  the  computational  capability  to  assess  future  exposure  concentra¬ 
tions,  enabling  efficient  analysis  of  tradeoffs  between  environmental  risk  (cleanup 
level)  and  costs  for  a  variety  of  cleanup  alternatives  prior  to  their  implementation. 

The  specific  technical  objectives  of  this  project  were  as  follows: 

a.  Develop  and/or  enhance  state-of-the-art  remediation  simulators  for  the 

following  technologies: 

( 1 )  Engineered  or  enhanced  in  situ  bioremediation. 

(2)  Natural  attenuation  of  petroleum  and  chlorinated  hydrocarbons  and 
explosives  (NAX). 

(3)  In  situ  chemical  treatment  (ISCT),  such  as  permeable  reactive  walls. 

(4)  Surfactant/cosolvent  flushing  to  recover  nonaqueous  phase  liquids 
(NAPLs). 

(5)  Venting  technologies,  such  as  soil  vapor  extraction  (SVE),  air 
sparging,  and  bioventing/sparging. 

b.  Validate  simulators  for  most  of  these  technologies  against  available 

laboratory  and  field  data. 

c.  Document  the  simulators  and  incorporate  them  into  the  GMS. 

The  purpose  of  this  report  is  to  summarize  the  results  achieved  during  the  course 
of  this  4-year  project. 


Technical  Approach 

Remediation  simulator  development  could  have  proceeded  along  three  paths, 
in  order  of  priority:  (a)  use  existing,  proven  remediation  simulators  where  avail¬ 
able  and  consistent  Avith  project  goals;  (b)  modify  promising  groundwater  model 
codes  to  simulate  additional  technologies  as  appropriate;  or  (c)  develop  new  codes 
as  required  for  efficient  simulation  of  innovative  technologies.  All  simulator  codes 
selected  or  developed  within  this  project  were  required  to  be  in  the  public  domain. 
Prior  to  this  project,  the  major  thrust  in  remediation  simulation  within  the  GMS 
had  been  toward  using  existing  codes  for  the  simpler  remediation  technologies, 
such  as  pump-and-treat  and  physical  barriers.  Codes  within  the  GMS  prior  to  this 
project  included  FEMWATER  (Lin  et  al.  1997),  MODFLOW  (McDonald  and 
Harbaugh  1988)  and  MT3D  (Zheng  1990).  The  approach  used  during  this  project 
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followed  (a)  and  (b),  i.e.,  to  use  existing  codes  and  to  modify,  adapt,  and  enhance 
existing  codes  to  simulate  additional  technologies  as  necessary. 

A  plan  for  GMS  integration  of  additional  simulators  existed  prior  to  this 
project  (Hadala  et  al.  1993)  that  placed  priority  on  certain  technologies  based  on 
remediation  technology  maturity  and  effectiveness,  user  demand  for  a  given 
remedial  technology,  and  the  degree  to  which  the  additional  simulator  would 
complement  existing  GMS  capabilities.  Based  on  these  guidelines,  a  suite  of 
three-dimensional,  groundwater/subsurface  simulation  codes  was  selected  as 
shown  in  Table  1.  All  of  these  codes  existed  prior  to  this  project,  but  were 
modified  during  this  project  to  better  represent  the  targeted  in  situ  remediation 
technologies. 


Table  1 

CU-1062  Simulators,  Technologies,  and  Validation  Sites 


Simulator 

Technology  Simulation 

Validated 

Validation  Site  || 

Sequential  Electron  Acceptor 
Model  in  3D  (SEAM3D) 

Engineered  in  situ 
bioremediation 

Natural  attenuation  of  petroleum 
hydrocarbons  (BTEX) 

Natural  Attenuation  Test 

Site  (NATS),  Columbus 

AFB 

Natural  attenuation  of 
chlorinated  hydrocarbons 
(solvents) 

Pensacola  Naval  Air 

Station 

University  of  Texas  Chemical 
Flood  Simulator  (UTCHEM) 

NAX 

Louisiana  Army 

Ammunition  Plant  (LAAP) 
in  collaboration  with  CU- 
715 

Surfactant  flushing 

Hill  AFB 

Cosolvent  flushing 

Lab  studies  and  Dover  AFB 

Operator  Splitting  in  3D 
(OS3D) 

ISCT,  permeable  reactive  walls 

Dover  AFB,  in  collaboration 
with  CU-107  and  Moffett 

Field 

Nonisothermal,  Unsaturated/ 
Saturated  Flow  and  Transport 

In  3D  (NUFT3D) 

SVE 

Biosparging 

Lawrence  Livermore 

National  Laboratory 

Building  518,  DOE  Hill  AFB 

.  . —  .  . 

Note:  SEAM3D  is  a  reaction  module  that  must  be  coupled  with  a  transport  module.  SEAM3D  is 
coupled  to  the  MT3DMS  transport  code,  which  resides  within  the  GMS  and  has  the  option  of  using 
either  SEAM3D  or  RT3D  for  reactions. 


The  four  targeted  simxilators  can  simulate  multiple  technologies.  For  example, 
UTCHEM  (University  of  Texas  Chemical  Flood  Simulator)  is  applicable  to  cosol¬ 
vent  flushing  and  surfactant  flushing.  Through  a  separate  EPA  effort,  UTCHEM 
has  been  modified  to  include  bioremediation  (aerobic  and  anaerobic  respiration 
and  cometabolism).  UTCHEM  is  a  multiphase  flow/transport  code,  so  it  poten¬ 
tially  offers  other  simulation  capabilities.  NUFT  (Nonisothermal  Unsaturated/ 
Saturated  Flow  and  Transport)  can  simulate  SVE,  steam-  or  fliermally-enhanced 
SVE,  air  sparging,  as  well  as  other  multiphase  flow/transport  problems.  NUFT  has 
also  been  modified  for  bioremediation.  There  is  some  overlap  of  simulator  capa¬ 
bilities,  such  as  multiphase  flow/transport  with  UTCHEM  and  NUFT  and  bio¬ 
remediation  simulation  with  SEAM3D  (Sequential  Electron  Acceptor  Model, 
3-Dimensional),  NUFT,  and  UTCHEM.  However,  these  overlaps  are  a  strength 
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rather  than  a  weakness  since  each  of  these  codes  has  its  unique  capabilities. 
Additionally,  model  selection  for  a  particular  application  should  be  based  on 
satisfying  study  needs  while  using  the  code  that  is  easiest  to  use. 

All  simulators  were  validated  against  available  laboratory  and  field  data,  with 
preference  given  to  Strategic  Enviromnental  Research  and  Development  Program 
(SERDP)  National  Environmental  Technology  Test  Sites  (NETTS).  Table  1  lists 
the  technology  simulations  that  were  validated  and  the  validation  sites.  Results  of 
all  validation  applications  either  have  been  or  are  in  the  process  of  being  docu¬ 
mented  as  technical  reports  or  papers,  with  references  cited  herein.  Each  simulator 
has  been  documented  and  implemented  in  the  GMS.  Also,  a  tutorial  for  each 
model  is  available  with  the  GMS  for  assistance  in  application. 

Before  any  model  is  used  for  decision  making,  it  should  be  validated  by  appli¬ 
cation  to  at  least  one  field  study.  Validation  provides  the  assurance  that  the  algo¬ 
rithms  embodied  in  the  computer  code  correctly  represent  the  physical  processes 
or  system  to  which  the  model  is  applied  (Beljin  1988).  A  model  is  said  to  be  vali¬ 
dated  when  sufficient  testing  shows  an  acceptable  degree  of  matching  the  actual 
systems  through  model  versus  data  comparisons.  However,  a  model  should  be 
revalidated  every  time  it  is  modified  or  is  applied  to  a  new  system.  In  practice, 
validation  is  a  never-ending  process  since  models  of  any  value  tend  to  stay  in  use 
for  a  long  time,  and  are  eventually  applied  to  a  wide  variation  of  field  systems. 
Additionally,  model  developers/users  are  always  striving  for  perfection  in  terms 
of  adding  new  features  to  make  models  more  closely  resemble  the  real  system. 
Therefore,  numerical  models  are  assessed  (validated)  by  repeated  application  to 
laboratory  experiments  and  monitored  field  data.  Repeated  successful  applications 
enhance  confidence  in  the  model,  while  unsatisfactory  results  stimulate  either 
improvements  to  the  model  or  its  abandonment  by  the  user  community.  While  the 
validation  applications  imdertaken  during  this  project  may  not  be  the  first  for  a 
particular  model,  they  most  likely  will  not  and  should  not  be  the  last. 

This  report  provides  a  summary  of  the  work  conducted  under  this  project.  The 
next  foxu"  chapters  focus  on  each  of  the  four  simulators  that  were  under  develop¬ 
ment.  Each  chapter  describes  the  simulator,  its  enhancements  or  improvements,  its 
validation,  its  incorporation  within  the  GMS,  and  recommendations  for  future 
work.  For  more  detail,  readers  should  refer  to  other  reports  and  papers  cited  in 
each  chapter. 
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Chapter  1  Introduction 


Figure  1 .  SEAM3D  incorporates  MT3DMS  packages  (blue)  for  advective- 

dispersive  transport,  sorption,  and  fluid  source/sinks  and  additional 
packages  unique  to  SEAM3D  (red),  including  the  Biodegradation, 
NAPL  Dissolution,  Cometabolism,  and  Reductive  Dechlorination 
Packages. 
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2  SEAM3D  and  MT3DMS 


Model  Description 

SEAM3D  is  a  computer  code  for  the  numerical  simulation  of  solute  transport 
with  aerobic  and  sequential  anaerobic  biodegradation.  The  model  can  depict 
multiple  constituents  in  a  three-dimensional,  anisotropic,  heterogeneous  domain. 
SEAM3D  is  designed  for  application  to  engineered  (accelerated)  bioremediation 
systems  and  intrinsic  bioremediation  (natural  attenuation). 

The  parent  code  is  MT3DMS  (Zheng  and  Wang  1 999),  which  is  designed  to 
simulate  physical  transport,  sorption,  and  mass  injection/extraction  due  to  fluid 
sources  and  sinks.  As  shown  in  Figure  1,  SEAM3D  consists  of  four  additional 
packages:  Biodegradation,  NAPE  Dissolution,  Cometabolism,  and  Reductive 
Dechlorination.  The  Biodegradation  and  NAPE  Dissolution  Packages  were 
developed  and  verified  in  the  original  version  of  SEAM3D  (Waddill  and 
Widdowson  2000).  The  latter  two  packages  were  recently  developed,  and  that 
work  is  outlined  in  this  report. 

The  Biodegradation  Package  is  the  main  subroutine  of  SEAM3D.  Hydro¬ 
carbon  contaminants  are  simulated  as  electron  donors  (i.e.,  substrates)  for 
microbial  growth,  with  available  electron  acceptors  (EAs)  used  in  the  following 
sequence  based  on  energy  yield:  oxygen  (O2),  nitrate  (NOp,  oxidized  manganese 

(Mn(IV)),  ferric  iron  (Fe(III)),  sulfate  (SO^” ),  and  carbon  dioxide  (CO2).  Mineral 

nutrients  and  biodegradation  end  products  and  daughter  products  are  also  simu¬ 
lated.  Biodegradation  of  each  substrate  follows  Monod  kinetics,  modified  to 
include  the  effects  of  EA  and  nutrient  availability.  Inhibition  functions  allow  any 
EA  to  inhibit  utilization  of  all  other  EAs  that  provide  less  energy  to  the  microbes. 
Microbial  biomass  is  simulated  as  scattered  microcolonies  attached  to  the  porous 
medium.  Equations  used  in  the  Biodegradation  Package  are  described  in  Waddill 
and  Widdowson  (1998). 

The  NAPE  Dissolution  Package  simulates  the  transfer  of  mass  from  a  NAPE 
to  the  aqueous  phase.  When  using  the  Biodegradation  Package  with  the  NAPE 
Dissolution  Package,  the  user  specifies  the  number  of  tracers  and  hydrocarbon 
substrates  and  the  characteristics  of  each  constituent,  and  the  characteristics, 
distribution,  and  concentration  of  a  light  nonaqueous  phase  liquid  (ENAPE)  such 
as  gasoline.  Equations  used  in  the  NAPE  Dissolution  Package  are  described  in 
Waddill  and  Widdowson  (1998). 


Chapter  2  SEAM3D  and  MT3DMS 


5 


Model  Enhancements 

Model  enhancements  to  SEAM3D  are  summarized  here.  Enhancements  to 
MT3DMS  can  be  found  in  Zheng  and  Wang  (1999). 


Reductive  dechlorination  package 

The  SEAM3D  Reductive  Dechlorination  Package  (SEAM3D-RDP)  is  a 
numerical  model  designed  to  simulate  intrinsic  bioremediation  of  chlorinated 
ethenes  in  groundwater.  SEAM3D-RDP  is  designed  to  simulate  transport  and 
biodegradation  of  chlorinated  compounds  perchloroethene  (PCE),  trichloroethene 
(TCE),  1,2-dichloroethene  (cis-DCE),  and  vinyl  chloride  (VC).  End  products 
simulated  in  SEAM3D-RDP  include  ethene  and  chloride.  The  model  is  solved 
numerically  using  the  core  packages  of  MT3DMS  and  the  Biodegradation 
Package.  SEAM3D-RDP  may  be  executed  with  or  without  the  NAPE  Dissolution 
and  Cometabolism  Packages.  The  NAPE  Dissolution  Package  was  updated  to 
simulate  PCE  or  TCE  dissolution  from  a  dense  nonaqueous  phase  liquid 
(DNAPE),  such  as  a  solvent.  When  the  Cometabolism  Package  and 
SEAM3D-RDP  are  used,  only  TCE,  cis-DCE,  and  VC  (or  any  combination) 
may  be  specified  as  subject  to  cometabolism. 

SEAM3D-RDP  is  designed  to  simulate  two  or  three  significant  biological 
degradation  mechanisms:  reductive  dehalogenation  (or  dechlorination)  and  direct 
oxidation.  A  third  mechanism,  cometabohsm,  is  simulated  m  the  Cometabolism 
Package,  if  this  option  is  specified  by  the  user.  The  mechanisms  for  degradation 
are  summarized  in  Table  2  for  the  four  chlorinated  ethenes.  The  conceptual  and 
mathematical  models  for  SEAM3D-RDP  and  relevant  equations  are  discussed  in 
depth  in  a  progress  report  to  ERDC  and  in  a  peer-reviewed  journal  paper  under 
preparation.  In  addition,  this  work  as  been  presented  in  part  at  two  international 
coifrerences  and  at  the  SERDP  2000  Symposium. 


Recent  experimental  studies  have  demonstrated  aerobic  and  anaerobic  oxida¬ 
tion  of  VC  and  cis-DCE,  but  direct  oxidation  appears  limited  to  higher  energy- 
yielding  terminal  electron  accepting  processes  (TEAPs).  All  four  simulated 


Table  2 

Biodegradation  Mechanisms  Considered  in  the  Model  for 
Perchloroethene  (PCE),  Trichloroethene  (TCE),  Dichloroethene 
(DCE),  and  Vinyl  Chloride  (VC) 

Compound 

Mechanism  I 

Direct  Oxidation 

Reductive 

Dechlorination 

Cometabolism 

PCE 

No 

Yes 

No 

TCE 

No 

Yes 

Yes' 

DCE 

Yes 

Yes 

Yes’ 

VC 

Yes 

Yes 

Yes’ 

^  If  specified  by  user  in  the  Cometabolism  Package. 
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chlorinated  ethenes  may  imdergo  reductive  dechlorination  with  the  rate  of  trans¬ 
formation  dependent  on  the  prevailing  TEAP  simulated  in  any  given  model  cell. 
SEAM3D  simulates  the  prevailing  TEAP  in  each  model  cell  and  at  every  time- 
step  using  the  Biodegradation  Package,  and  this  information  is  linked  to  the 
Reductive  Dechlorination  Package.  VC  and  cis-DCE  may  serve  as  electron  donors 
in  the  SEAM3D-RDP  when  conditions  in  a  model  cell  favor  direct  oxidation. 


Cometabolism  package 

The  SEAM3D  Cometabolism  Package  (SEAM3D-CP)  is  a  numerical  model 
designed  to  simulate  engineered  or  intrinsic  bioremediation  of  recalcitrant  com¬ 
pounds  in  groundwater.  SEAM3D-CP  is  designed  to  simulate  transport  and 
biodegradation  of  chlorinated  ethenes  (TCE,  cis-DCE,  trans-DCE,  VC)  or  other 
compoxmds  known  to  cometabolize  in  the  presence  of  a  single  growth  substrate, 
oxygen,  and  an  aerobic  microbial  population.  The  model  is  solved  nmnerically 
using  the  core  packages  of  MT3DMS  and  the  Biodegradation  Package. 
SEAM3D-CP  may  be  used  in  conjunction  with  the  NAPE  Dissolution  Package, 
SEAM3D-RDP,orboth. 

When  SEAM3D-CP  is  run,  one  hydrocarbon  compound  is  designated  in  the 
Biodegradation  Package  as  the  growth  substrate  linked  to  cometabolism.  A  frac¬ 
tion  of  the  aerobic  population  (designated  in  the  Biodegradation  Package)  is 
assumed  responsible  for  cometabolism.  Methane,  an  end  product  of  aerobic  bio¬ 
degradation,  may  also  serve  as  the  growth  substrate.  In  this  case,  a  methanotrophic 
population,  separate  from  the  main  aerobic  population,  is  designated  in  the 
Biodegradation  Package.  Up  to  three  recalcitrant  compounds  may  be  designated  in 
the  SEAM3D-CP.  Equations  used  in  the  SEAM3D-CP  were  modified  from 
Semprini  and  McCarty  (1991  and  1992). 


Microbial  populations 

The  SEAM3D  Biodegradation  Package  originally  included  as  many  as  six 
bacterial  populations:  aerobes,  facultative  nitrate  reducers,  anaerobic  manganese 
reducers,  anaerobic  iron  reducers,  anaerobic  sulfate  reducers,  and  methanogens 
(Table  3).  Three  new  microbial  groups  have  been  added.  Two  biomass  concen¬ 
tration  variables  are  simulated  in  SEAM3D-RDP,  both  responsible  for  the  reduc¬ 
tion  of  chlorinated  ethenes.  One  population  utilizes  only  PCE  and  TCE  as  EAs, 
preferentially  using  PCE,  but  inhibited  by  higher  energy-yielding  EAs.  The 
second  population  utilizes  only  cis-DCE  and  VC  as  EAs,  preferentially  using  cis- 
DCE,  but  potentially  inhibited  in  the  presence  of  PCE  and  TCE.  As  noted  in  the 
preceding  paragraph,  a  seventh  microbial  population,  the  methanotrophs,  was 
added  directly  to  the  Biodegradation  Package.  This  group  is  relevant  only  to  the 
SEAM3D-CP,  when  cometabolism  is  simulated.  The  methanotrophic  population 
may  be  simulated  either  with  or  without  the  methanogenic  population. 
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Table  3 

Electron  Acceptors  Used  by  the  Microbial  Populations  for  Biodegradation  of 
carbon  Substrates 

fdro- 

Microbial  Population 
x  =  SEAM3D  Bio 
Y=SEAM3D’RDP 

Index  of  Electron 

Donor  Is 

Index  of 
Electron 
Acceptor  /e^ 

End 

Products 

x=^  strict  aerobes 

User-specified  hydrocarbon 
(UHC) 

/e=1  O2 

- 

x=2  facultative  NO3  reducers 

UHC 

/e=1  O2 

- 

UHC 

/e=2  NO3 

O2 

NO, 

x=3  anaerobic  Mn(IV)  reducers 

UHC 

/e=3  Mn(x) 

O2,  NO3,  Mn(IV),  Fe(lll),  SO4 

Mn(li) 

x=4  anaerobic  Fe(lll)  reducers 

UHC 

/e=4  Fe(x) 

O2,  NO3,  Mn(l\/),  Fe(lll),  SO4 

Fe(ll) 

x=5  anaerobic  SO4  reducers 

UHC 

/e=5  SO4 

O2,  NO3,  Mn(IV),  Fe(lll),  SO4 

HS- 

x=6  methanogens 

UHC 

O2,  NO3,  Mn(IV),  Fe(lll),  SO4 

CH4 

x=7  methanotrophs 

Methane 

/e=1  O2 

- 

- 

y=1  PCE/TCE  reducers 

Not  simulated 

O2,  NO3,  Mn(IV),  Fe(lll),  SO4 

Cl 

y=2  cis-DCEA/C  reducers 

Not  simulated 

O2,  NO3,  Mn(IV),  Fe(lll),  SO4 

Cl,  ethene 

’  Electron  acceptors  are  listed  In  order  of  highest  to  lowest  Gibbs  free  energy  per  half  reaction. 

2  Utilization  of  each  electron  acceptor  is  Inhibited  by  the  presence  of  an  electron  acceptor  that  provides  higher  energy. 

Model  Validation 

Waddill  and  Widdowson  (2000)  demonstrated  the  use  of  SEAM3D  at  a 
gasoline-contaminated  site  where  the  uncontrolled  release  resulted  in  a  benzene, 
toluene,  ethylbenzene,  xylene  (BTEX)  plume  that  included  methyl  tertiary  butyl 
ether  (MTBE).  The  model  was  validated  using  monitoring  well  data  collected 
over  5  years.  The  SEAM3D  Biodegradation  and  NAPE  Dissolution  Packages 
were  more  recently  validated  at  the  Natural  Attenuation  Study  (NATS)  located  at 
Columbus  Air  Force  Base  (AFB),  Mississippi.  Progress  was  also  made  toward 
validation  of  the  SEAM3D-RDP  at  the  Pensacola  Naval  Air  Station  (NAS). 


NATS 

NATS  consisted  of  the  release  and  subsequent  monitoring  of  a  NAPE  mixture 
of  BTEX,  naphthalene,  decane,  and  bromide  into  the  shallow,  unconfined  aquifer 
underlying  Columbus  AFB  in  eastern  Mississippi.  During  the  present  study,  both 
conceptual  and  mathematical  models  were  developed  for  NAPE  source  release, 
sequential  aerobic/anaerobic  biodegradation,  and  sorption  during  NATS. 
SEAM3D  was  used  to  simulate  fully  three-dimensional  transport  and  aerobic, 
anoxic,  ferrogenic,  and  methanogenic  hydrocarbon  biodegradation.  Simulation 
results  matched  individual  BTEX  concentration  distributions  collected  5  and 
9  months  following  NAPE  release. 

NATS  was  initiated  by  placing  a  known  mass  of  NAPE  into  the  Columbus 
aquifer.  The  mass  firaction  of  BTEX  and  naphthalene  in  the  NAPE  was  propor¬ 
tional  to  a  typical  IP-4  jet  fuel,  with  decane  serving  as  an  octane  surrogate  for  the 
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remaining  inert  fraction  of  the  fuel.  The  NAPL  was  mixed  with  30  m^  of  clean 
sand  and  placed  in  an  excavated  trench  (hereafter  referred  to  as  the  source  trench), 
which  was  backfilled  to  land  surface  with  imcontaminated  sand.  Shee^iling, 
which  had  been  installed  along  the  walls  of  the  trench  to  establish  hydraulic 
control,  was  removed,  and  the  natural  groundwater  gradient  in  the  source  area  was 
restored  to  allow  petroleum  hydrocarbon  (PHC)  dissolution  into  the  groruidwater. 
Solute  concentrations  were  measured  from  selected  ports  in  a  multilevel  sampling 
network.  The  most  comprehensive  collection  and  analysis  were  conducted  during 
the  sampling  events  5  months  and  9  months  after  source  removal.  Analysis  of 
solute  concentrations  indicated  that  multiple  microbial  processes  were  active 
during  the  NATS  experiment. 

Figure  2  presents  contour  plots  of  field  data  (left  column)  at  time  t  =  5  months 
and  calibrated  single  porosity  model  results  (right  column)  for  a  vertical  cross 
section  located  4  m  west  of  the  source  center  line  for  benzene,  toluene,  ethyl¬ 
benzene,  p-xylene,  and  naphthalene.  Comparison  between  columns  in  Figure  2 
shows  that  the  calibrated  model  recreated  the  fingering  observed  for  the  1-mg/L 
BTEX  contour,  and  generally  reproduced  the  vertical  and  longitudinal  extent  of 
BTEX  concentration  plumes.  The  vertical  and  longitudinal  dimensions  of  the 
naphthalene  plume  were  also  reproduced,  but  field-measured  naphthalene  concen¬ 
trations  in  and  around  the  source  area  were  approximately  five  times  higher  than 
model-predicted  values  at  t  =  5  months  (possibly  due  to  increased  solubility 
effects).  A  comparison  of  model  results  att  =  9  months  (results  not  shown) 
showed  a  stable  plume  and  the  persistent  fingering  on  the  hydrocarbon  plume. 

Horizontal  slices  of  the  field-measured  and  model-simulated  PHC  concen¬ 
trations,  taken  at  an  elevation  z  that  passed  through  the  source  trench  (i.e.,  60  m), 
showed  that  the  reactive  transport  model  also  captured  horizontal  variations  in 
PHC  transport.  As  an  example.  Figure  3  depicts  horizontal  contour  plots  of 
toluene  concentration  developed  from  field  data  (left  column)  and  the  calibrated 
single  porosity  model  (right  column)  for  i  =  5  months  and  9  months.  Figure  3 
shows  agreement  between  field  observations  and  model  simulations  in  that 
toluene  migration  in  the  western  portion  (negative  transverse  distance)  of  the 
plume  was  greater  than  in  the  eastern  portion. 

SEAM3D  mass  balance  calculations  indicated  that  biodegradation  consumed 
49  percent  of  the  aqueous  phase  hydrocarbon  compounds,  while  13  percent  of  this 
mass  was  sorbed  and  38  percent  remained  in  the  aqueous  phase.  Mass  calculations 
further  indicated  that  aerobic  biodegradation  accounted  for  the  majority  of  hydro¬ 
carbon  attenuation  (37  percent  of  the  attenuated  mass),  followed  by  ferrogenesis 
(23  percent),  sorption  (21  percent),  nitrate  reduction  (16  percent),  and  methan- 
ogenesis  (<  4  percent).  Model  results  were  particularly  sensitive  to  the  NAPL 
release  rate,  the  initial  Fe(III)  concentration,  hydrocarbon  utilization  rates,  initial 
condition  for  the  anaerobic  microbial  populations,  and  dispersivity.  Results  of  this 
study  have  helped  to  confirm  the  validation  of  SEAM3D  for  simulating  the  trans¬ 
port  and  attenuation  of  petroleum  hydrocarbon  compounds  in  groundwater.  The 
results  of  the  SEAM3D  validation  for  the  NATS  are  documented  in  Braimer  and 
Widdowson  (accepted  for  publication).  Validation  of  MT3DMS  for  the  NATS  is 
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Figure  2.  Comparison  of  field-measured  BTEX  and  naphthalene  concentrations  with  model  simulations  for  the  single  porosity  model  along  a 
vertical  cross  section  4  m  west  of  the  source  center  line  for  f  =  5  months 


Figure  3.  Comparison  of  field-measured  toluene  concentrations  (left)  with  model-simulated  toluene 
concentrations  (right)  for  a  horizontal  cross  section  at  z  =  60  m  at  f  =  5  months  (top)  and 
9  months  (bottom) 

documented  in  Julian  et  al.  (in  review).  Additionally,  the  transport  scheme  in 
MT3DMS  is  documented  and  verified  in  Zheng,  Wang,  and  Dortch  (in  review). 


NAS  Pensacola 

Model  validation  was  confirmed  at  NAS  Pensacola  Wastewater  Treatment 
Plant  (WWTP),  where  natural  attenuation  of  chlorinated  ethene  compounds  has 
been  monitored  and  documented.  Past  disposal  of  industrial  wastewater  into  an 
unlined  holding  pond  led  to  contamination  of  the  surficial  aquifer  underlying  the 
WWTP.  Trichloroethene  is  present  with  a  concentration  range  of  2.5  mg/L  to 
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nondetect  along  the  91 -m  (300-ft)  flow  path.  Chlorinated  daughter  products  are 
also  present,  including  cis-DCE  and  VC.  Groundwater  monitoring  data  demon¬ 
strate  the  removal  of  all  three  chlorinated  ethenes  in  wells  immediately 
up-gradient  of  Pensacola  Bay. 

Marine  and  fluvial  terrace  sediments  underlie  the  WWTP,  consisting  pri¬ 
marily  of  fine  to  medium  sands  to  a  depth  of  12  m  (40  ft)  below  land  surface. 

The  unconfined  aquifer  is  underlain  by  a  6-m-  (20-ft-)  thick  confining  layer.  An 
upward  hydraulic  gradient  prevents  downward  transport.  Local  groundwater  flow 
is  to  the  west,  which  results  in  discharge  to  the  adjacent  Pensacola  Bay,  91  m 
(300  ft)  to  the  west  of  the  source  area.  Geochemical  data  indicate  that  iron- 
reducing  and  sulfate-reducing  conditions  are  the  primary  redox  processes  in  the 
aquifer  but  that  methanogenesis  is  also  active. 

The  modeling  objective  was  to  calibrate  the  model  to  contaminant  concen¬ 
tration  data  (TCE,  cis-DCE,  and  VC)  and  geochemical  parameters  along  the  main 
transect  of  the  plume.  A  flow  and  transport  model  of  the  NAS  Pensacola  WWTP 
site  was  developed  to  simulate  the  transport  and  attenuation  of  TCE,  cis-DCE,  and 
VC  using  MODFLOW  and  SEAM3D,  respectively.  Reduction  in  the  source  con¬ 
centration  of  TCE  due  to  chemical  oxidation  treatment  is  being  simulated  using 
the  calibrated  site  model  as  a  verification  step. 

Figures  4  and  5  show  TCE  and  cis-DCE  with  VC  concentration,  respectively, 
versus  distance  along  the  plume  center  line  for  both  the  measured  and  simulated 
cases.  Results  of  the  calibrated  model  show  a  stable  TCE  plume  with  an  accurate 
representation  of  the  down-gradient  peaks  of  cis-DCE  and  VC.  The  site  model 
using  SEAM3D  will  be  used  to  simulate  the  residual  TCE  plume  following  source 
remediation  using  chemical  oxidation.  SEAM3D  will  be  used  to  determine  the 
time  required  for  the  TCE  plume  to  restabilize  following  chemical  oxidation  of 
the  source  and  the  time  to  reach  remediation  goals  under  a  monitored  natural 
attenuation  program.  The  final  results  of  the  NAS  Pensacola  validation  will  be 
documented  in  a  paper. 


Interface  with  GMS 

A  revised  SEAM3D  interface  to  GMS  has  been  developed  that  will  incor¬ 
porate  the  Reductive  Dechlorination  and  Cometabolism  Packages.  In  addition,  the 
“define  species”  menu  has  been  revised  to  incorporate  chlorinated  ethenes  and 
methanotrophs.  The  NAPL  Dissolution  Package  input  was  also  revised  to  account 
for  solvents  derived  from  a  DNAPL. 


Recommendations 

Recommendations  for  further  validation  work  and  improvements  to  the 
SEAM3D  code  include  the  following: 
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Figure  4.  Measured  and  simulated  {SEAM3D-RDP)  concentration  of  TCE  along  the  main  plume 
transect  showing  a  steady-state  distribution  for  each  solute  prior  to  chemical  oxidation 
treatment  of  source  area  (to  convert  distance  into  meters,  multiply  by  0.3048) 

a.  Validate  the  SEAM3D  Cometabolism  Package  using  field  data  firom  a 
pilot  test  or  full-scale  bioremediation  project. 

b.  Validate  the  SEAM3D  Reductive  Dechlorination  Package  at  a  field  site 
where  a  biostimulant  (e.g.,  vegetable  oil  or  formate)  has  been  injected  in 
the  source  area  to  enhance  reductive  dechlorination  of  PCE/TCE  in  the 
source  area. 

c.  Develop  and  test  a  SEAM3D  Bioaugmentation  Package  and  validate  at  a 
controlled  pilot  test  or  full-scale  field  demonstration. 

d.  Develop  and  test  a  SEAM3D  Package  for  modeling  with  inhibition 
kinetics,  which  could  be  validated  at  the  Laurel  Bay,  SC,  U.S.  Geological 
Survey  study  site. 
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Figure  5.  Measured  and  simulated  (SEAM3D-RDP)  concentration  of  cis-DCE  and  VC  along  the  main 
plume  transect  showing  a  steady-state  distribution  for  each  solute  prior  to  chemical  oxidation 
treatment  of  source  area  (to  convert  distance  into  meters,  multiply  by  0.3048) 
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3  UTCHEM 


Model  Description 

UTCHEM  is  a  three-dimensional,  finite  difference,  multiphase  flow,  multi- 
component,  reactive  transport  model.  This  is  a  mature,  multipurpose  model  with 
too  many  facets  to  detail  here.  Major  features  are  summarized,  followed  by 
discussions  of  the  enhancements  made  during  the  current  project,  benchmarking 
and  validation  of  the  model,  development  of  the  GMS  interface,  and  the  prepara¬ 
tion  of  documentation  and  training  materials. 

The  key  UTCHEM  features  include  the  following: 

a.  Multidimensional:  fully  three-,  two-,  or  one-dimensional  or 
nondimensional  (batch)  in  Cartesian,  radial,  or  stretched  vertical 
coordinate  systems. 

b.  Multiphase  flow:  water,  oil,  air,  microemulsion;  the  aqueous  phase  is 
either  water  or  microemulsion  depending  on  local  surfactant 
concentration. 

c.  Multicomponent  reactive  transport:  any  number  of  user-defined,  organic 
and/or  inorganic  contaminants,  microbes,  or  additives  (e.g.,  reactants, 
nutrients,  etc.). 

d.  Constitutive  models  for  relative  permeability-saturation-capillary 
pressure  (kr-S-Pc)  relations:  Brooks-Corey  Imbibition  and  First-drainage 
models. 

e.  Hysteretic  models  (Parker-Lenhard)  for  capillary  pressure  and  relative 
permeability  in  water-wet  media  and  for  two-phase  flow  (oil-water)  in 
mixed-wet  media. 

/  Liquid  phase  densities  functions  of  composition. 

g.  Liquid  phase  viscocities  functions  of  composition,  including  viscosity¬ 
enhancing  polymers  with  non-Newtonian  rheology. 

h.  Full  tensor  coefficients  for  dispersivity  and  molecular  diffusion  in  each 
phase. 

i.  Microemulsion  phase  behavior:  empirical  models  dependent  on  complex 
function  of  surfactant/cosolvent,  oil  composition,  water  chemistry. 
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j.  Surfactant/cosolvent  effects  on  constitutive  models,  i.e.,  enhanced  oil 
mobilization  (versus  dissolution)  as  a  function  of  viscous  and  buoyancy 
forces  (trapping  number). 

k.  Partitioning  interwell  tracer  tests  (FfTTs):  characterization  of  oil  satura¬ 
tions  in  the  vadose  or  phreatic  zones  by  assessing  the  relative  transport  of 
alcohols,  varying  in  their  hydrophobicity  (preferential  partitioning  into 
any  oil  phase). 

/.  Wells:  flexible  module  for  any  number  of  injection,  withdrawal,  or 
monitor  wells,  completed  (or  screened)  in  user-specified  cells,  in  any 
orthogonal  orientation. 

m.  Equilibrium  or  rate-limited  mass  transfer  (dissolution)  of  NAPL 
components. 

n.  Equilibrium  adsorption  to  aquifer  solids  (linear,  Langmuir,  or  Freundlich 
isotherms);  cation  exchange. 

o.  Boundary  conditions:  impermeable  by  default;  user  may  specify  head  and 
compositional  conditions  at  two  opposing  boundaries;  flow  may  be 
confined  or  imconfined. 

p.  Initial  conditions:  user  may  specify  initial  phase  saturations,  solute 
concentrations. 

g.  Zero-order,  first-order  (a  common,  empirical  approach),  and  higher-order 
biogeochemical  models  accommodated. 

r.  Bioremediation  package:  highly  flexible  in  terms  of  reaction  pathways 
and  kinetics  (first-order  and  hi^er),  microbial  biomass  (microcolony 
concept). 

s.  Geochemistry  package:  limited  accommodation  of  user-defined, 
equilibrium  reactions,  e.g.,  speciation,  precipitation-dissolution. 

t.  Heterogeneous  media:  parameters  for  flow  and  transport  models  can  be 
assigned  as  homogeneous,  stratified,  or  fully  heterogeneous  fields. 

u.  Dual  porosity  approximation  of  fractured  media. 

V.  Heat  transport  (polythermal  conditions)  coupled  to  other  properties  (e.g., 
viscosity). 

w.  Implicit  pressme  explicit  saturation  (IMPES)  type  formulation;  third- 
order-accurate  finite-difference  methods  with  total-variation-diminishing 
(TVD)  flux  limiter. 

Detailed  documentation  and  user  information  regarding  all  of  these  features  are 
available  (Reservoir  Engineering  Research  Program  2000a,  2000b).  The  new 
GMS  interface  will  expedite  the  preparation  of  input  for  most  UTCHEM  9.0 
features,  particularly  those  relevant  to  environmental  applications.  The  remaining 
features,  for  which  a  GMS  interface  is  not  developed,  still  can  be  used,  but  will 
require  editing  of  the  input  file  independent  of  the  GMS. 

UTCHEM  employs  a  three-dimensional,  block-centered,  finite-difference 
scheme.  The  solution  method  is  implicit  in  pressure  and  explicit  in  concentration 
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(an  IMPES  type  method).  Options  for  one-  and  two-point  upstream  and  third- 
order  spatial  tessellation  are  available.  A  TVD  flux  limiter  has  been  added  to 
increase  the  stability  and  robustness  (reducing  overshoot  and  undershoot)  of  the 
second-  and  third-order  methods.  Of  course,  the  third-order  method  gives  the  most 
acciuate  solution. 

UTCHEM  was  developed  originally  in  the  late  1970’s  (Pope  and  Nelson 
1978)  as  a  design  simulator  for  surfactant-polymer  flooding  in  enhanced  oil 
recovery.  The  model  continued  to  evolve  and  improve  in  terms  of  both  the 
processes  described  and  numerical  methods  used.  One  of  the  essential  require¬ 
ments  of  a  smfactant  enhanced  aquifer  remediation  (SEAR)  design  simulator  is 
accurate  and  flexible  treatment  of  microemulsion  phase  behavior,  which  is  central 
to  the  treatment  of  other  phenomena  such  as  partitioning,  interfacial  tension, 
capillary  pressure,  capillary  number,  and  microemulsion  viscosity.  Adaptation  of 
the  mature  and  successful  modeling  approach  to  enhanced  oil  recovery  evident  in 
UTCHEM  was  deemed  a  rational  and  cost-effective  approach  to  the  development 
of  a  SEAR  simulator.  Some  modifications  were  required  in  the  adaptation  to 
environmental  concerns,  arising  primarily  from  the  focus  on  much  lower  oil 
saturations  (subresidual)  and  very  low  concentrations  of  organic  contaminants  in 
the  aqueous  phase  (parts  per  billion  (ppb)  level).  These  adaptations  have  been 
affected  under  efforts  supported  recently  by  the  EPA  and  the  current  SERDP 
Project  CU-1062.  Two  goals  of  the  current  project  were  to  develop  a  graphical 
user  interface  within  the  GMS  and  to  continue  the  evolution  of  UTCHEM  as  an 
advanced  remediation  simulator. 

UTCHEM  is  a  flexible,  multifaceted  model  that  would  be  useful  in  the  design 
and  evaluation  of  a  variety  of  remediation  technologies.  However,  the  current 
effort  is  focused  on  smfactant  and  cosolvent  flushing  for  NAPE  cleanup.  The 
validation  process  was  accomplished  by  assessment  of  the  revised  model  in  appli¬ 
cation  to  several  surfactant/cosolvent  flushing  experiments  at  the  laboratory  and 
field  scales.  Equally  important  was  the  thorough  beta  testing  of  the  GMS  interface 
as  it  developed.  Several  of  these  applications  will  serve  as  tutorials  in  the  self- 
directed  training  of  new  UTCHEM-GMS  users. 


Model  Enhancements 

The  enhancements  to  UTCHEM  under  SERDP  CU-1062  have  been  achieved 
in  three  basic  areas:  (a)  expanding  model  capabilities  for  environmental  applica¬ 
tions,  particularly  for  SEAR-related  processes,  (b)  developing  a  graphical  user 
interface  to  improve  ease  of  use,  and  (c)  improving  numerical  methods.  Each  of 
these  enhancement  areas  is  discussed  more  thoroughly  in  this  section.  Applica¬ 
tions  demonstrating  these  enhancements  are  presented  in  the  “Model  Validation” 
section. 

Several  enhancements  to  UTCHEM  are  either  “invisible”  to  the  user  or 
involve  processes  that  could  be  useful  in  the  design/evaluation  of  multiple  remedi¬ 
ation  technologies.  The  former  would  include  the  shift  to  dynamic  memory  alloca¬ 
tion  and  other  changes  in  numerical  methods  to  take  advantage  of  FORTRAN  90 
features.  The  dynamic  memory  allocation  allows  array  sizes  to  adjust 
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automatically  to  match  the  input,  thereby  avoiding  the  need  to  redimension  the 
arrays  explicitly  and  recompile  the  entire  model  for  each  application.  Enhance¬ 
ments  that  could  benefit  the  design  and  evaluation  of  multiple  remediation 
technologies  include  the  following: 

a.  The  addition  of  Freimdlich  and  Langmuir  adsorption  isotherms. 

b.  Adaptation  of  input  and  output  file  formats  for  GMS  compatibility. 

c.  Expansion  of  boimdary  condition  options. 

Independent  efforts,  concurrent  with  CU-1062,  also  have  contributed  to 
UTCHEM  enhancement,  which  are  incorporated  in  UTCHEM-GMS  (ver¬ 
sion  9.0).  These  enhancements  supported  by  DOE,  the  petroleum  industry,  or 
others  corroborate  the  broad  appeal  and  utility  of  this  simulator. 

The  numerous  enhancements  to  UTCHEM  necessitated  an  update  to  the 
documentation  (Reservoir  Engineering  Research  Program  2000a,  2000b).  The 
revised  user’s  guide  includes  descriptions  of  new  input  parameters  and  formats. 
The  technical  documentation,  which  is  a  thorough  suite  of  brief  reviews  of 
imderlying  theory  and  mathematical  descriptions  of  key  UTCHEM  features,  is 
revised  to  reflect  the  more  fimdamental  changes.  Guidance  is  offered  as  to  param¬ 
eter  determination  and  typical  values  and  ranges  for  select  parameters,  with  an 
emphasis  on  those  unique  to  SEAR.  The  revised  user’s  guide  and  techmcal  docu¬ 
mentation  accompany  any  release  of  the  GMS  interface  and  currently  are  available 
in  PDF  format  from  the  UTCHEM  page  at  the  GMS  Web  site, 
http://chl.wes.armv.mil/software/gms/docs.htD. 

Validation  of  the  model  and  testing  of  the  GMS  interface  will  be  documented 
in  separate  reports  and/or  integrated  into  the  GMS  documentation  as  tutorials.  The 
reports  will  include  discussion  of  the  rationale  for  parameter  selection,  particularly 
for  parameters  unique  to  surfactant/cosolvent  flush  simulation.  Finally,  the  GMS 
documentation  for  UTCHEM,  to  accompany  each  public  release,  will  incorporate 
portions  of  each  of  these  documents,  with  an  emphasis  on  the  graphical  interface 
features. 


Enhancements  to  surfactant/cosolvent  flushing  simulation 

UTCHEM  is  already  well  established  as  one  of  the  most  comprehensive 
models  available  for  SEAR  design  and  evaluation.  In  addition  to  the  GMS  inter¬ 
face,  there  are  several  SEAR-related  enhancements  including  the  following: 

a.  Nonequilibrium  mass  transfer  to  the  micellar  phase. 

b.  Hysteretic  constitutive  models  (relative  permeability-capillary  pressure). 

c.  Constitutive  models  for  mixed-wet  media  (for  oil-water  systems). 

All  of  these  enhancements  were  completed  and  tested  in  Fiscal  Year  (FY)  1999. 
The  GMS  interface  will  serve  only  to  accelerate  the  acceptance  of  both  UTCHEM 
and  smfactant/cosolvent-flushing  technologies  by  raising  the  level  of  under¬ 
standing  in  the  user  community  and  providing  the  requisite  tools  to  design  and 
implement  SEAR.  UTCHEM  also  provides  a  tool  for  the  forward  modeling 
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assessment  of  PITTs  for  estimating  NAPL  saturations  prior  to  SEAR  and  in  the 
subsequent  assessment  of  its  effectiveness. 

Most  experienced  SEAR  practitioners  would  concur  that  simulation  is  critical 
in  the  design  and  preinstallation  evaluation  phase  and  is  not  a  peripheral  or 
academic  exercise.  At  the  very  least,  design/evaluation  simulation  is  the  best 
approach  for  the  quantitative  integration  of  multiple,  coupled,  physical  and  chem¬ 
ical  processes  that  can  affect  SEAR  performance.  There  is  no  alternative  for  the 
integration  of  complex  calculations  of  compositional  effects  on  such  factors  as 
phase  equilihria  or  constitutive  models  for  flow  in  a  heterogeneous  medium. 


Enhancements  to  the  bioremediation  features  of  UTCHEM 

Natural  or  engineered  microbial  processes  can  be  described  or  designed  using 
the  highly  flexible,  input-controlled  definition  of  reaction  pathways  and  mecha¬ 
nisms  available  in  UTCHEM.  These  features  may  be  useful  in  the  evaluation  of 
the  post-SEAR  cleanup  of  any  residual  NAPL  or  contaminated  groundwater,  as 
well  as  any  unrecovered  surfactant/cosolvent. 

The  monitored  natural  attenuation  of  organic  contaminants  can  be  an  attrac¬ 
tive,  cost-effective  cleanup  alternative  where  tenable.  The  attenuation  of  certain 
recalcitrant  contaminants,  such  as  solvents  and  explosives,  can  involve  slow, 
subtle  processes  imder  the  predominantly  aerobic,  oligotrophic  conditions  present 
in  most  aquifers.  Thus,  long-term  fate  and  transport  predictions  become  indispen¬ 
sable  in  the  evaluation  of  monitored  natural  attenuation  efficacy.  Groundwater 
modeling  provides  the  only  rigorous  approach  to  consider  and  evaluate  monitored 
natural  attenuation  alternatives  and  assess  exposure  risk. 

UTCHEM  has  been  made  more  readily  applicable  to  the  simulation  of  natural 
attenuation  and  the  design  of  engineered  bioremediation.  The  major  enhance¬ 
ments  to  UTCHEM  bioreaction  features  under  CU-1062  included  adding  the 
following: 

a.  A  first-order  kinetic  model  alternative  that  accommodates  product 
generation. 

b.  A  microcolony-type,  bioreaction  kinetics  model. 

c.  A  built-in  transport  time-step  control,  based  on  dimensionless  parameters 
conditioned  upon  the  relative  influence  of  reaction  kinetics  versus 
diffusive  mass  transfer  on  bioreaction  rates. 

d.  Option  to  consider  the  coupled  effect  of  microbial  biofouling  (growth)  on 
flow. 

e.  Option  to  specify  initial  conditions  for  aqueous  phase  transport. 

These  features  build  on  the  bioreaction  module  developed  for  EPA  (UTCHEM 
version  6.0). 

Bioreaction  kinetics  can  be  influenced  by  solute  mass  transfer  between  pore 
water  and  biomass.  The  error  associated  with  neglecting  mass  transfer  limitations 
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to  biomass  is  a  function  of  both  the  rate  of  biodegradation  and  the  rate  of  mass 
transfer.  The  rate  of  biodegradation  also  influences  the  operator  splitting  error  (a 
small  error  inherent  in  the  solution  scheme),  so  that  it  is  difficult  to  distinguish 
between  errors  introduced  from  neglecting  mass  transfer  and  errors  introduced  by 
operator  splitting.  To  resolve  this  difficulty,  the  following  tasks  were  completed: 
(a)  the  dependence  of  the  operator  splitting  error  on  the  flow  rate  and  biodegrada¬ 
tion  rate  has  been  established  in  terms  of  dimensionless  numbers,  and  confirmed 
with  a  large  number  of  one-dimensional  simulations;  and  (b)  UTCHEM  will  auto¬ 
matically  select  the  appropriate  biodegradation  kinetics  time-step,  independent  of 
the  advection-dispersion  time-step.  The  Damkohler  and  Peclet  dimensionless 
numbers  are  used  to  determine  when  diffusive  mass  transfer  to/from  biomass  is 
significant.  Generally,  a  Damkohler  number  below  0.1  indicates  that  mass  transfer 
is  not  limiting  the  bioreactions.  The  automatic  time-step  feature  greatly  reduces 
simulation  times  by  solving  the  biodegradation  equations  only  as  often  as  neces¬ 
sary  to  keep  the  operator  splitting  error  beneath  a  given  threshold.  This  approach 
also  allows  better  comparisons  between  runs  with  and  without  mass  transfer. 

Biofouling  of  porous  media  by  overstimulating  microbial  growth  is  a  phenom¬ 
enon  to  be  avoided  in  bioremediation  design.  The  biofouling  algorithm  incorpo¬ 
rated  into  UTCHEM  is  based  on  porosity  loss  due  to  a  constrained  biomass 
growth.  The  altered  permeability  is  expressed  using  a  Kozeny-Carmen  equation. 
Biomass  growth  is  limited  to  occupying  up  to  90  percent  of  the  available  pore 
space.  Biomass  transport  and  adsorption  also  have  been  added  to  UTCHEM  with 
the  adsorption  modeled  as  linear  partitioning. 

UTCHEM  has  been  applied  to  a  variety  of  bioreaction  problems,  including 
evaluation  of  explosives  reactions  at  the  Louisiana  Army  Ammunition  Plant  and 
bioremediation  of  post-SEAR,  residual  solvent  and  surfactant  at  the  Hill  AFB 
pilot  demonstration  site  (de  Blanc  1998).  One  or  more  of  these  will  serve  as 
tutorials  for  the  GMS  documentation. 


Model  Validation 

A  variety  of  conditions  have  been  simulated  with  UTCHEM  9.0  to  validate 
the  model  and  to  evaluate/test  its  GMS  interface.  The  immediate  goal  was  to  ferret 
out  as  early  as  possible  the  coding  errors  that  are  inevitable  in  models  as  complex 
as  UTCHEM  and  GMS.  Validation  of  the  more  unique  and  advanced  process 
models  within  UTCHEM,  most  of  which  are  not  amenable  to  exact  (analytical) 
solution,  requires  model  application  to  a  variety  of  well-controlled  and/or 
well-monitored  experiments  at  laboratory  and  field  scales.  Repeated  success 
enhances  confidence  in  the  model,  while  unsatisfactory  results  stimulate  either 
improvements  to  the  model  or  its  abandonment  by  the  user  community.  This 
natural  selection  phenomenon  in  model  evolution  is  driven  by  Are  same  market 
pressures  at  work  in  the  remediation  business  in  which  the  better,  faster,  cheaper 
technologies  survive. 


‘  C.  J.  McGrath  and  M.  Zakikhani.  (1998).  “Numerical  modeling  of  the  natural  attenuation  of 
explosives,”  presentation  at  3rd  Tri-Service  Environmental  Technology  Workshop,  San  Diego,  CA, 
18-20  August  1998. 
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UTCHEM  has  been  demonstrated,  benchmarked,  and  validated  in  application 
to  a  variety  of  scenarios,  including  (a)  surfactant  flushing  for  the  removal  of 
residual  DNAPL  at  Hill  AFB,  Utah;  (b)  PITTs  at  Hill  AFB;  (c)  ethanol  flushing  to 
recover  residual  benzene  from  packed  columns,  (d)  an  alcohol  flushing  demon¬ 
stration  for  DNAPL  removal  at  Dover  AFB,  and  (e)  post-SEAR  bioremediation  at 
Hill  AFB.  Each  of  these  UTCHEM  applications  was  supported,  at  least  in  part, 
under  SERDP  CU-1062,  while  the  experimental/field  data  were  collected  for 
other  projects. 


Surfactant  flushing,  Hill  AFB 

UTCEDEM  simulations  were  indispensable  in  the  design  of  the  highly  success¬ 
ful  surfactant  flush  demonstration  at  Hill  AFB,  Utah  (Brown  et  al.  1996).  These 
simulations  have  been  revised  to  reflect  the  flow  conditions  actually  employed 
(versus  designed)  and  with  an  updated  format  for  UTCHEM  and  GMS  (Martino, 
in  preparation).  The  purpose  was  to  provide  a  well-documented,  three- 
dimensional  example  of  UTCHEM  application  to  SEAR  to  be  included  with  the 
code  documentation.  This  history  matching  exercise  generally  improved  the  agree¬ 
ment  between  predicted  and  observed  data.  Most  model  parameters,  such  as  those 
related  to  phase  properties,  constitutive  models,  hydrogeology,  and  the  numerical 
grid  were  unchanged. 

The  PITT  represents  an  innovative  method  to  estimate  in  situ  NAPL  satura¬ 
tions  and  possibly  the  distribution  (Figure  6).  Saturations  are  estimated  based  on 
the  relative  retardation  of  tracers  of  contrasting  partitioning  coefficient  Kp,  with 
more  hydrophilic  tracers  being  effectively  conservative  (e.g.,  isopropyl  alcohol 
(IP A),  Kp  -  0),  while  relatively  hydrophobic  tracers  (e.g.,  1-heptanol,  Kp  =  140.5; 
pentanol,  Kp  =  3.9)  are  more  severely  retarded  in  the  presence  of  NAPL.  PITTs 
were  used  to  assess  the  initial  DNAPL  saturations  at  the  Hill  site  (Jin  1995;  Jin  et 
al.  1997;  Duke  Engineering  and  Services  1998).  A  second  (“intermediate”)  PITT 
was  conducted  after  a  brief  initial  surfactant  flush  (Phase  I).  A  third  (“final”)  PITT 
was  performed  after  a  second,  more  substantial  surfactant  flush  (Phase  II).  The 
similarity  in  arrival  times  for  IPA,  pentanol,  and  heptanol  in  the  final  PITT 
suggests  the  absence  of  significant  NAPL  saturation. 

UTCHEM  can  simulate  the  effects  of  adding  water-soluble  polymers  or 
forming  foam  to  further  enhance  the  effectiveness  of  SEAR.  The  addition  of  a 
water-soluble,  food-grade  polymer  to  enhance  viscosity  or  air  injection  to  induce 
foaming  increases  the  sweep  efficiency  of  the  flushing  solution.  Greater  sweep 
efficiency  enhances  the  recovery  of  NAPL  that  otherwise  might  be  bypassed  and 
left  behind  in  lenses  of  lower  permeability.  This  capability  of  UTCHEM  has  been 
demonstrated  for  creosote  removal  (Wu  et  al.  2000). 


Ethanol  flushing 

Alcohol  flushing  for  the  recovery  of  residual  NAPL  from  porous  media  has 
been  documented  in  the  laboratory  and  under  controlled  field  demonstrations. 
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Figure  6.  UTCHEM  simulations  of  PITT  and  SEAR  at  Hill  AFB  OU-2: 

(A)  Predicted  versus  observed  PITT  for  2-propanol  breakthrough  at 
well  SB-1;  (B)  Predicted  versus  observed  PITT  for  1-pentanol;  note 
that  breakthrough  of  the  more  lipophilic  pentanol  {Kp  =  3.9)  is  retarded 
slightly  relative  to  IPA  {Kp  =  0);  and  (C)  Simulated  oil  saturations 
before  and  after  SEAR 
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Though  the  procedure  is  technically  feasible,  issues  of  efficiency  and  cost  remain 
to  be  addressed.  Lower  molecular  weight  alcohols  such  as  methanol  and  ethanol 
serve  to  enhance  the  solubility  of  NAPL  contaminants  in  aqueous  solution,  thus 
the  term  cosolvent. 

UTCHEM  has  been  applied  to  a  series  of  column  experiments  designed  to 
elucidate  critical  flow  and  transport  phenomena  involved  in  ethanol  flushing  for 
residual  NAPL  (benzene)  recovery  (McGrath,  in  preparation).  Several  concentra¬ 
tions  of  ethanol  in  aqueous  solution  were  flushed  through  a  40-cm  column  of  veiy 
fine  sand  (~0. 10  ±  0.05  mm)  containing  residual  benzene  saturation  (residual 
saturation  Sro  =  0.20).  At  low  ethanol  concentration  (0.50  volume  fraction  (VF)), 
benzene  dissolution  was  enhanced,  but  no  mobilization  was  evident.  At  0.67  VF 
ethanol,  benzene  solubility  into  the  aqueous  phase  was  enhanced  as  expected  and 
a  small  amount  of  benzene  was  mobilized  as  an  unstable  macroemulsion  (Fig¬ 
ure  7).  At  0.75  VF  ethanol,  benzene  mobilization  was  sufficiently  effective  to 
form  a  distinct  oil  bank  in  advance  of  the  ethanol  front.  Magnetic  resonance 
imaging  (MRI),  used  to  monitor  the  advance  of  the  fluids  through  the  columns, 
revealed  some  interesting  displacement  phenomena  (Figure  7b).  The  effluent  was 
analyzed  by  pulse  nuclear  magnetic  resonance  (NMR). 

Simulation  of  the  alcohol  flush  experiments  with  UTCHEM  (Figure  7a,  c) 
reproduced  most  of  the  flow  phenomena  that  this  model  could  be  reasonably 
expected  to  capture,  such  as  enhanced  dissolution,  NAPL  mobilization,  the 
buoyant  override  of  the  alcohol.  Effluent  flow  rates  and  bulk  concentrations  are 
reproduced  quite  well.  However,  observed  phenomena  such  as  viscous  fingering 
and  the  formation  of  macroemulsions  are  random  or  nonthermodynamic  phenom¬ 
ena  that  are  too  poorly  understood  to  predict/model  accurately.  The  use  of 
UTCHEM  for  alcohol  flushing  requires  a  few  ad  hoc  redefinitions  of  model 
parameters  that  are  intended  for  surfactant  flushing.  For  example,  the  requisite 
definition  of  the  critical  micelle  concentration  (CMC)  has  no  meaning  in  alcohol 
flushing.  If  the  alcohol  is  modeled  essentially  as  a  poor  surfactant,  as  has  been 
done  here,  a  “pseudo-CMC”  may  be  redefined  by  the  user  as  the  alcohol  concen¬ 
tration  above  which  the  aqueous  solubility  of  NAPL  constituents  is  substantially 
enhanced.  Such  compromises  are  not  unusual  when  pushing  the  limits  of  a  still 
evolving  code  such  as  UTCHEM.  These  compromises  are  not  problematic  as  long 
as  the  limitations  are  understood.  However,  a  more  general  cosolvent  model  that 
captures  the  full  range  of  concentration-dependent  partitioning  could  be  devel¬ 
oped  for  UTCHEM  or  any  alcohol  flush  simulator.  Additionally,  UTCFDEM 
assumes  negligible  cosolvent  (surfactant)  solubility  in  the  oil  phase  (i.e.,  plait 
point  is  set  at  the  oil  apex),  which  is  generally  not  the  case  with  low  molecular 
weight  alcohols  and  definitely  not  the  case  for  larger,  more  hydrophobic  alcohols. 
Again,  this  appears  to  be  a  minor  limitation,  but  a  more  general  liquid-liquid 
equilibria  model  would  be  preferable. 

The  UTCHEM  simulator  continues  to  be  evaluated  further  in  application  to 
the  cosolvent  flush  demonstration  recently  completed  at  Dover  AFB.  The  Dover 
demonstration  involves  the  controlled  release  of  a  known  amount  of  DNAPL 
(Figure  8)  within  a  laterally  contained  section  of  natural,  mildly  heterogeneous, 
sandy  aquifer.  An  extensive,  high-resolution  (spatially  and  temporally). 
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Figure  7.  Cosolvent  flush  (0.67  VF  ethanol)  of  NAPL  (benzene)  removal  from  a 
horizontal  packed  column:  (A)  UTCHEM-predicted  oil  saturation 
distribution  at  0.5  pore  volume  (PV)  flush.  (B)  MRI  of  column  section  at 
approximately  0.5  PV  flush;  note  the  near-vertical  advanced  front, 
followed  by  an  axial,  convex-forward  front,  similar  to  the  predicted 
phenomena.  (C)  UTCHEM-predicted  (dashed  lines)  versus  observed 
(points)  effluent  flow  rate  (upper  plot)  and  total  composition  deter¬ 
mined  by  pulse  NMR 
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Figure  8.  UTCHEM  simulation  of  PCE  release  prior  to  the  cosolvent  flush  demonstration  at  the  Dover 
AFB  test  cell.  A  total  of  92  L  of  PCE  was  introduced  via  12  wells  (eight  9-L  and  four  5-L 
injections)  within  the  4.6-  by  3.0-m  test  cell 


three-dimensional  data  set  for  both  preflush  and  postflush  PITTs  as  well  as  the 
alcohol  flush  are  providing  the  benchmark  against  which  to  evaluate  UTCHEM 
as  an  alcohol  flush  simulator. 


Post-SEAR  bioremediation  at  Hill  AFB 

The  bioremediation  features  of  UTCHEM  have  been  demonstrated  in  the 
design  simulation  of  cometabolic  bioremediation  of  postsurfactant  flush  residuals 
after  the  SEAR  demonstration  at  Site  0-2,  Hill  AFB,  Utah  (de  Blanc  1998).  The 
intent  was  to  demonstrate  these  new  model  features  and  verify  that  they  are  fully 
functional  and  predict  reasonable  results.  Injection  of  dissolved  methane  and 
oxygen  (peroxide)  was  simulated  to  assess  the  cometabolic  biodegradation  of  the 
TCE  by  indigenous  methanotrophic  microorganisms.  Initial  quantities  of  TCE 
included  approximately  5  L  of  residual  TCE  (of  the  ~1,890  L  present  initially) 
and  dissolved  TCE  leaching  into  the  SEAR-treated  zone  from  adjacent  untreated 
areas.  Dissolved  oxygen  in  the  aquifer  was  depleted  initially  due  to  heterotrophic 
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degradation  of  residual  surfactant.  The  TCE  is  presumed  to  have  suppressed  the 
indigenous  methanotroph  population.  Three  injection  wells  and  three  recovery 
wells  were  used  to  establish  a  plug  flow.  Pump-and-treat  was  considered  as  a 
baseline  scenario  for  which  oxygen-saturated  water  was  injected  for  200  days. 
Much  of  the  surfactant  was  removed,  but  TCE  remained  high  due  to  continued 
dissolution  of  remaining  NAPE  traces  and/or  diffusion  from  low-permeability 
zones.  Injection  of  methane-saturated  (~20  mg/L),  highly  oxygenated  (300  mg/L 
as  hydrogen  peroxide)  water  rejuvenated  the  methanotrophs,  particularly  near  the 
injection  well,  where  TCE  concentrations  were  reduced  due  to  the  cometabolic 
biodegradation.  Unfortunately,  the  methane  was  consumed  quickly,  precluding 
significant  penetration  into  the  aquifer,  compared  to  a  nonreactive  tracer.  The 
surfactant  concentrations  were  also  lower  than  those  for  the  pump-and-treat 
scenario  due  to  the  higher  oxygen  influx.  Although  cost  savings  over  pump-and- 
treat  could  be  realized  by  this  methane-peroxide  approach,  performance  might  be 
enhanced  still  further  by  considering  alternative  scenarios  such  as  groundwater 
reinjection  and  cycled  injection  of  methane  and  peroxide.  This  simple  example 
demonstrates  one  of  the  advantages  of  modeling  during  the  design  phase  to  screen 
viable  alternatives  for  the  most  promising. 


GMS  Implementation 

A  three-phase  approach  was  adopted  for  the  development  of  the  GMS  inter¬ 
face  to  UTCHEM,  compatible  with  the  comprehensive  but  modular  structure  of 
the  model.  Phase  1,  completed  early  in  FY  99,  yielded  an  interface  to  the  grid 
module  only.  The  Phase  2  release  (FY  00)  included  an  interface  to  the  essential 
input/output  variables  for  surfactant  and/or  cosolvent  flushing.  The  third  and  final 
phase  (to  be  released  in  2001)  includes  development  of  the  interface  to  the 
partitioning  interwell  tracer  features  and  the  comprehensive  and  flexible  bio¬ 
remediation  package. 

The  first  UTCHEM-GMS  Workshop  was  held  on  16-18  November  1998  at 
the  University  of  Texas  at  Austin.  This  workshop  provided  an  opportunity  for  the 
interface  developers,  UTCHEM  users  (new  and  old),  and  SEAR  practitioners  to 
exchange  ideas  regarding  the  model  and  interface.  Over  25  groundwater  profes¬ 
sionals  from  DoD,  DOE,  EPA,  academia,  and  private  consulting  firms  were 
invited  to  participate.  Both  the  numerical  aspects  of  UTCHEM  and  the  biogeo¬ 
chemical  processes  involved  in  SEAR  and  other  cleanup-related  features  of 
UTCHEM  were  discussed.  Several  model  applications  were  presented,  including 
SEAR,  alcohol  flushing,  and  natural  attenuation  of  explosives,  followed  by  a  half¬ 
day  GMS  training  session. 

All  releases  of  UTCHEM-GMS  will  be  accompanied  by  the  tutorials 
described  in  this  section  and  additional  tutorials  as  they  become  available.  Tuto¬ 
rials  expedite  the  learning  of  a  new  code  by  example  and  are  particularly  useful 
when  learning  comprehensive  simulators  such  as  UTCHEM.  These  tutorials  are 
designed  to  lead  the  UTCHEM  and/or  GMS  novice  through  applications,  intro¬ 
ducing  additional  features  with  each  tutorial. 
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Tutorial-1  (T-01) 


The  first  tutorial  simulates  a  PCE  spill  in  a  simple,  structured  heterogeneous 
media  (Figure  9).  This  tutorial  reproduces  the  familiar  DNAPL  experiments  of 
Kueper  (1989)  familiar  to  most  investigators  of  DNAPL  flow  phenomena.  This 
tutorial  demonstrates  the  following: 

a.  Setting  up  a  uniform  numerical  grid. 

b.  Assigning  heterogeneous  media  properties  for  multiphase  flow. 

c.  Definition  of  a  simple  GMS  map  overlay  (the  media  interfaces). 

The  tutorial  also  demonstrates  the  numerical  stability  of  UTCFEM.  Figure  9 
shows  the  PCE  saturation  predicted  by  UTCHEM  at  310  sec  after  the  initiation  of 
PCE  release. 


Tutorial-2  (T-02) 


The  second  tutorial  simulates  a  PCE  spill  and  redistribution  followed  by  a 
long  cleanup  by  pump-and-treat  (Figure  10).  A  steady  release  of  28.32  m^ 
(1,000  ft^)  of  PCE  near  the  surface  of  a  horizontally  stratified,  sandy  medium  is 
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Figures.  UTCHEM  simulation  of  DNAPL  (PCE)  flowthrough  heterogeneous  media  for  Tutorial-1.  Media 
boundaries  are  outlined  by  dashed  lines.  PCE  was  injected  at  the  center  top  and  allowed  to 
flow  for  310  sec  (5.17  min).  This  GMS  tutorial  based  on  experiments  by  Kueper  (1989) 
demonstrates  the  stability  of  UTCHEM  in  heterogeneous  media 
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A.  30-days  of  PCE  release 


G.  Day  90  H.  Day  96 

(SPAR  day  0)  (SEAR  day  6) 


I.  Day  100 
(SEAR  day  10) 


J.  Day  102  K.  Day  104 
(SEAR  day  22)  (SEAR  day  24) 


Figure  10.  UTCHEM  simulations  for  UTCHEM-GMS  Tutorial-2  (water  flush)  and 
Tutorial-3  (surfactant  flush),  based  upon  the  Borden  SEAR  demon¬ 
stration.  Images  are  DNAPL  saturations  (PCE)  using  a  uniform  scale. 
(A)  Distribution  after  30  days  of  PCE  release  (0.943  m^/day 
(33.3  ft^/day));  the  central  withdrawal  well  is  flanked  by  two  injection 
wells;  regional  flow  gradient  is  left  to  right.  (B  and  G)  Distribution  after 
60  days  of  redistribution.  (B-F)  Time  series  of  water  flush,  analogous 
to  a  pump-and-treat  scenario;  PCE  is  dissolved  after  ~1 5  years.  (G-K) 
Time  series  of  surfactant  flush;  removal  is  completed  in  approximately 
25  days.  Note  that  surfactant  from  the  left  well  reaches  the  recovery 
well  on  day  94,  while  surfactant  from  the  right  well  (further  away  and 
down  gradient)  does  not  break  through  until  day  102,  thus  the 
asymmetry  and  delay  beyond  10  days  of  surfactant  flush 


simulated.  The  PCE  is  allowed  to  redistribute  for  60  days  within  a  moderate 
regional  gradient,  reaching  the  aquifer  base  to  form  a  small  pool.  DNAPL 
saturations  are  up  to  about  0.27  at  day  90.  A  pump-and-treat  scheme  is  imposed 
using  one  fully  penetrating  recovery  well  beneath  the  spill  location,  flanked  by 
two  injection  wells  (at  about  9  and  16  m  (30  and  52  ft)  away).  The  trapped 
DNAPL  dissolves  slowly,  being  most  persistent  in  the  four  lower-permeability 
horizons  and  a  basal  pool.  The  PCE  saturations  reach  zero  by  day  5600,  i.e., 

15  years  of  pump-and-treat. 
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Tutorial-3  (T-03) 


The  third  tutorial  builds  upon  T-02  by  replacing  the  lengthy  pump-and-treat 
period  with  a  16-day  surfactant  flush  to  achieve  complete  DNAPL  removal 
(Figure  10).  This  UTCHEM  simulation  introduces  the  user  to  the  additional  input 
parameters  required  for  SEAR  simulation.  Tutorials  T-02  and  T-03  are  based  on 
the  SEAR  demonstration  at  the  Borden  site  and  UTCHEM  simulations  by  Brown 
(1993)  and  Jin  (1995).  In  practice  this  scenario  would  probably  use  a  surfactant 
flush  much  less  than  the  16  days  simulated,  and  may  involve  a  pre- surfactant  slug 
of  high-salinity  water  to  condition  the  aquifer  for  optimal  SEAR  effectiveness. 

The  research  site  at  Canadian  Forces  Base,  Borden,  Ontario,  subsequently 
referred  to  as  the  Borden  site,  is  well  characterized  in  terms  of  aquifer  hydro¬ 
geology,  solute  transport  properties,  and  multiphase  flow  properties.  Several  tracer 
dispersion  and  remediation  technology  demonstration  studies  have  been  con¬ 
ducted  at  the  Borden  site,  including  the  first  demonstration  of  zero-valent  iron  for 
in  situ  treatment  of  dissolved  chlorinated  organic  solvents.  The  thin  surficial/ 
unconfmed  aquifer  consists  of  relatively  homogeneous,  medium-  to  fine-grained, 
siliciclastic  sands  underlain  by  a  thick  clay  unit  at  a  depth  of  ~9  m  (~30  ft).  The 
water  table  is  2  to  3  m  (6.5  to  10  ft)  below  the  ground  surface.  The  hydraulic 
conductivity  is  —0.00723  cm/sec  (—20.5  ft/day),  typical  for  clean  sands. 

Additional  tutorials  in  preparation  will  highlight  many  of  the  other  useful 
features  of  UTCHEM  and  the  GMS  interface.  Specifically,  these  tutorials  include 
demonstration  of  the  following:  (a)  partitioning  interwell  tracers  with  which  to 
characterize  NAPE  saturations  in  either  saturated  or  unsaturated  media;  (b)  the 
bioremediation  features  and  how  they  may  be  coupled  to  the  fate  of  NAPE  con¬ 
taminants;  and  (c)  cosolvent  flushing  for  NAPE  removal.  Still  more  tutorials  will 
be  prepared  in  the  future  as  the  capabilities  of  UTCHEM  continue  to  expand. 


Recommendations 

UTCHEM-GMS  brings  a  powerful  and  practical  remediation  simulator  to  the 
GMS  suite.  UTCHEM-GMS  is  an  advanced,  multifaceted  model  well  suited  to 
the  following  common,  environmentally  significant  groundwater  problems: 

a.  NAPE  (dense  or  light)  spills  into  saturated  or  unsaturated  media. 

b.  NAPE  cleanup  by  flushing  with  surfactant/cosolvent/polymer  or 
surfactanl/foam. 

c.  PITTs  in  saturated  or  unsaturated  media. 

d  Bioremediation  by  natural  or  engineered  microbial  processes. 

e.  Geochemical  reactions  (e.g.,  heavy  metals,  radionuclides). 

There  are  very  few  multiphase-flow,  multicomponent  transport  models  in  the 
public  domain,  and  still  fewer  are  as  well  documented  and  maintained  as 
UTCHEM.  The  capabilities  for  simulating  NAPE  flow,  multicomponent  transport, 
PITTs,  SEAR,  and  bioremediation  within  a  single  model  make  UTCFfEM-GMS 
unique  among  remediation  simulators.  The  goal  of  the  GMS  graphical  interface  is 
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to  reduce  any  impedance  to  users  talcing  advantage  of  this  advanced  remediation 
simulator. 

UTCHEM,  like  all  models,  has  its  limitations  and  idiosyncrasies,  rooted  in  its 
developmental  history.  UTCHEM  was  developed  initially  to  describe  microemul¬ 
sion  phase  behavior  in  petroleum  reservoirs.  The  terminology  and  umts  used  in 
UTCHEM  reflect  this  heritage.  Though  grormdwater  modelers  without  a  chemical 
or  petroleum  engineering  background  may  find  this  distracting,  the  inconvenience 
is  minor. 

UTCHEM  was  not  originally  developed  for  cosolvent-only  flushing,  general 
vapor  phase  transport,  geochemical  reactions,  or  vadose  zone  problems;  but  due  to 
a  series  of  modifications  over  the  past  several  years,  the  code  can  be  applied  to  a 
wide  variety  of  conditions.  However,  for  some  types  of  applications,  users  may 
find  it  necessary  to  use  some  approximations.  For  example,  cosolvent  behavior 
can  be  approximated  as  a  poor  surfactant.  Limitations  in  application  to  cosolvent- 
only  flushing  will  be  addressed  in  a  separate  report.  Vapor  phase  transport  can  be 
simulated  for  the  special  case  of  air-oil  partitioning  tracers;  thus,  UTCHEM 
(version  9)  should  not  be  applied  to  the  simulation  of  vapor  extraction,  as  the 
NUFT  code  can.  A  geochemical  package  has  been  added  to  UTCHEM;  however, 
this  package  is  not  presently  part  of  the  GMS  interface.  UTCHEM  should  be 
applicable  to  vadose  zone  problems,  although  this  feature  was  not  validated 
during  this  project. 

Models  evolve  as  the  user  community  develops  a  more  complete,  quantitative 
understanding  of  the  processes  being  modeled  and/or  the  need  arises  to  extend 
model  capabilities  to  describe  new  processes.  Continuous  development  to  address 
new  model  needs  is  the  mark  of  a  vibrant  or  “living”  code;  failure  to  evolve  leads 
inevitably  to  abandonment  of  an  outdated  model.  UTCHEM  likely  will  continue 
to  be  improved  to  take  advantage  of  new  ideas. 

The  GMS  interface  developed  under  the  current  project  does  not  provide 
graphical  access  to  all  UTCHEM  capabilities.  The  neglected  features  are  those 
either  designed  for  petroleum  reservoir  applications  or  those  deemed  to  be  of 
lesser  utility  in  environmental  applications.  These  features  remain  in  the  model 
(version  9.0),  but  utilization  would  require  editing  of  the  input  files  outside  the 
GMS.  Specialized  UTCHEM  features  such  as  horizontal  wells,  radial  grids,  and 
stretch  grid  options  are  also  not  included  in  the  current  GMS  interface. 

The  geochemical  reaction  package  is  the  only  UTCHEM  package  entirely 
missing  fi-om  the  current  GMS  interface  (version  3.1).  The  geochemical  package 
is  somewhat  limited  and  geared  toward  chromium  cross-linking  reactions  impor¬ 
tant  in  surfactant-polymer  floods  for  enhanced  oil  recovery.  Another  GMS  model 
(OS3D)  is  designed  specifically  for  environmental  geochemical  reactive  transport 
modeling,  e.g.,  in  situ  chemical  treatment. 

The  current  versions  of  UTCHEM  and  GMS  visualize  the  distribution  of  fluid 
phase  saturations  and  total  concentration  for  each  component.  It  may  also  be  use¬ 
ful  in  the  design  phase  to  be  able  to  visualize  the  distribution  of  other  parameters 
and  conditions  such  as  phase  flow  velocities  or  relative  permeabilities,  total 
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component  concentrations  in  each  phase  separately,  interfacial  tension,  and 
trapping  number.  UTCHEM  calculates  and  outputs  all  of  the  variables  mentioned, 
but  the  output  is  not  presently  formatted  for  GMS  visualization.  Future  improve¬ 
ments  to  the  GMS  interface  for  UTCHEM  should  consider  adding  these  addi¬ 
tional  output  options. 

The  enhancements  to  UTCHEM  imder  CU-1062  advance  the  practice  of  both 
multiphase  flow  modeling  and  evaluation  and  design  of  SEAR  as  a  remediation 
technology.  UTCHEM  continues  to  be  an  excellent  vehicle  for  the  quantitative 
description  and  evaluation  of  remediation  processes  that  may  be  coupled  to 
synergistic  or  competing  processes.  GMS  interface  development  should  continue 
to  complete  graphical  input-output  to  the  ever-expanding  model  options  in 
UTCHEM. 
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Model  Description 

NUFT  is  a  general-purpose  computer  code  for  modeling  multiphase  fluid  flow 
and  multispecies  reactive  transport  in  porous  media  under  nonisothermal  or  iso¬ 
thermal  conditions  (Nitao  1998a).  It  has  been  used  for  various  field  applications  in 
enviromnental  remediation:  soil  vapor  extraction  (Nitao,  Martins,  and  Ridley 
2000),  soil  vapor  extraction  combined  with  groundwater  dewatering  (Rueth  et  al. 
1998),  bioventing  (Sim  et  al.  2000),  dynamic  steam  stripping  and  contaminant 
hydropyrolysis  (Newmark  et  al.  1999),  and  electrical  heating  (Camgan  and  Nitao 
2000).  NUFT  has  also  been  used  in  field  studies  for  research  in  vadose  zone  flow 
and  transport  processes  (Lee  and  Nitao  2000;  Carrigan  1999),  as  well  as  for  evalu¬ 
ation  of  vadose  zone  contaminated  sites  (Demir  et  al.  1999).  Other  field  applica¬ 
tions  include  nuclear  waste  disposal  (Nitao  and  Buscheck  1995),  nuclear  treaty 
verification  (Carrigan  et  al.  1996),  containment  of  gases  during  subcritical  test 
explosions,  and  enhanced  petroleum  recovery  (Sahni,  Kumar,  and  Knapp  2000). 

NUFT  is  a  robust  code  that  can  solve  highly  computationally  demanding 
problems.  It  runs  on  IBM  PC-compatibles  and  various  Unix  operating  systems 
such  as  Sun  Solaris,  IBM  ADC,  HP-UX,  DEC  OSF,  and  386  Linux.  The  input  file 
format  is  highly  flexible.  The  code  is  written  in  the  C++  language.  Because  it  uses 
dynamically  allocated  memory,  the  code  does  not  have  to  be  recompiled  for 
different  problem  dimensions. 

NUFT  consists  of  several  modules  contained  in  a  single  source  code.  Each 
module  has  its  own  set  of  simplifying  assumptions  so  that  the  user  can  select  the 
most  physically  appropriate  mathematical  module  and  computationally  efficient 
numerical  solution  method.  The  following  are  current  modules: 

a.  UCSAT  =  unconfined  aquifer  flow  model. 

b.  US  IP  =  Richard’s  equation  flow  model. 

c.  USIC  =  single-component,  single-phase  transport  model. 

d  USNT  =  general  multiphase,  multicomponent  flow  and  transport  under 
either  isothermal  or  nonisothermal  conditions. 

e.  JOULE  =  electrical  heating. 
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The  distinct  modules  in  the  code  use  a  common  set  of  utility  routines  and 
input  file  format.  Because  the  various  modules  are  essentially  isolated  fi'om  each 
other,  fiiture  modules  can  be  added  without  affecting  existing  ones.  Input  data  are 
in  the  form  of  that  used  by  the  LISP  language.  An  internal  LISP  interpreter  for  the 
Scheme  dialect  of  LISP  is  part  of  the  simulator  to  parse  and  store  input  data  infor¬ 
mation,  to  perform  input  error  checking,  and  to  set  default  input  parameter  values. 

USNT  is  the  NUFT  module  (Nitao  1998b)  that  is  the  focus  of  this  report.  The 
USNT  module  solves  the  multiphase  flow  and  multispecies  flow  and  transport 
equations  under  nonisothermal  or  isothermal  conditions.  Local  thermodynamic 
equilibrium  for  partitioning  of  species  between  fluid  phases  is  assumed.  Vapor 
pressme  lowering  of  components  is  available  as  an  option.  A  multiporosity  model 
is  also  available  as  an  option  for  modeling  fractured  porous  rock  systems.  Kinetic 
reaction  rate  laws  can  be  turned  on,  such  as  first-order,  power  law,  sequentially 
first-order,  Monod,  and  dual-substrate  Monod  reactions. 

Unlike  some  other  multiphase  codes,  the  NUFT-USNT  module  can  handle 
full  disappearance  of  the  NAPL,  or  any  fluid  phase,  due  to,  for  example,  dissolu¬ 
tion  or  evaporation.  This  is  in  contrast  to  some  codes  that  require  at  least  a  small 
amount  of  NAPL  to  be  present  everywhere  in  the  domain,  which  means  that,  for 
those  models,  the  predicted  groundwater  concentrations  can  never  go  below  the 
contaminant  solubility. 

The  NUFT-USNT  module  solves  the  partial  differential  equations  for  the 
conservation  of  mass  and  energy  using  the  integrated  finite-difference  method. 

The  resulting  nonlinear  system  of  equations  is  solved  at  each  time-step  by  the 
Newton-Raphson  method.  Each  iteration  of  the  Newton-Raphson  method  requires 
the  solution  of  a  system  of  linear  equations  arising  fi'om  the  linearization  of  the 
nonlinear  system.  Options  for  solution  of  the  linear  equation  system  are  the  vari¬ 
ably  banded  Gaussian  elimination  method  and  the  iterative  conjugate  gradient 
method  with  various  options  for  preconditioners.  See  the  NUFT  Reference 
Manual  (Nitao  1998a)  and  the  NUFT  USNT  User’s  Manual  (Nitao  1998b)  for  the 
actual  balance  equations  and  references  to  numerical  methods. 


Model  Enhancements 

Monod-type  reactions  were  implemented  as  part  of  the  kinetic  reaction 
capability  in  the  NUFT-USNT  module  to  model  bioremediation  of  contaminants, 
such  as  by  bioventing,  and  to  model  natural  degradation.  The  general  form  of  the 
implemented  reaction  law  is  a  product  of  factors  of  the  form  bJ(bi+Ci),  CJ(Si+C^), 
and  Cf'  where  C,  is  the  concentration  of  the  reactant  species,  h,  is  the  inhibition 
constant,  s,  is  the  saturation  constant,  and  p,  is  some  power.  The  index  i  can  run 
over  any  set  of  contaminants  or  other  mass  components  such  as  oxygen  or 
methane.  An  additional  temperature-dependent  factor  proportional  to  the 
Arrhenius-t5q)e  rate  law  exp(-F'//?Z)  is  also  optionally  present.  Here,  R  is  the  gas 
constant,  T  is  absolute  temperature,  and  E  is  activation  energy,  J/mol.  Documen¬ 
tation  of  the  Monod  rate  laws  using  NUFT  is  given  by  Sim  et  al.  (2000). 
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Another  enhancement  in  NUFT  was  the  capability  to  generate  a  hetero¬ 
geneous  field  of  saturated  and  unsaturated  hydrological  properties.  This  feature 
allows  for  more  realistic  modeling  of  remediation  sites  where  natural  heter¬ 
ogeneity  plays  an  important  role.  In  the  approach  described  in  this  report,  the 
physical  domain  is  subdivided  into  any  number  of  lithological  material  types.  Up 
to  three  statistically  independent  random  fields  are  generated  for  each  material 
type.  Hydrologic  properties  in  a  particular  material  are  specified  as  functions  of 
these  multiple  random  fields  so  that  the  generated  air  entry  pressure  field,  for 
example,  can  be  statistically  correlated  with  the  saturated  permeability  field.  The 
random  fields  are  generated  by  the  spectral  method  (Lee  and  Nitao  2000)  using 
the  fast  Fourier  transform  algorithm.  Figure  1 1  shows  an  example  of  random 
fields  generated  for  a  model  of  a  steam  flooding  test  performed  at  a  creosote- 
contaminated  site  at  Visalia,  CA  (Newmark  et  al.  1999).  Here  three  random  fields 
were  generated,  one  for  each  zone.  The  air  entry  pressure  parameter  for  the  capil¬ 
lary  pressure  curves  were  also  varied  spatially  by  assuming  that  they  scale  as  the 
square  root  of  the  permeability. 


Figure  1 1 .  Heterogeneous  random  fields  generated  for  saturated  intrinsic 
permeability  for  a  model  of  a  steam  flooding  test  performed  at  a 
creosote-contaminated  site  at  Visalia,  CA 

A  multicontinua  option  was  implemented  into  NUFT  that  enables  the  code, 
for  example,  to  model  sites  with  fractured  porous  rock  by  treating  the  fracture 
network  and  the  porous  rock  matrix  as  separate  porous  medium  continua.  Theo¬ 
retically,  any  number  of  continua  can  be  specified  to  model  separate  fracture  sets. 
Exchange  of  fluids  between  continua  is  governed  by  various  flow  options.  The 
random  field  generation  also  works  with  this  capability.  In  particular,  each  con¬ 
tinuum  can  have  material  types  with  each  material  type  having  up  to  three  statisti¬ 
cally  independent  random  fields.  The  file  addendum.doc  distributed  with  the 
NUFT  code  documents  the  multicontinua  and  random  property  field  options. 
Figure  12  shows  the  isosurface  of  50  percent  liquid  water  saturation  (saturation  is 
the  fractional  occupied  void  space)  in  the  fractures  due  to  infiltrating  water  within 
a  heterogeneous  domain. 
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Figure  12.  Water  infiltrating  within  a  heterogeneous  fracture  continuum  where  the 
50  percent  liquid  water  saturation  isosurface  is  shown 


Model  Validation 

A  primary  objective  of  this  project  was  to  validate  NUFT  for  SVE.  Code 
enhancements  and  general  support  provided  by  this  project  contributed  to  the 
preparation  and  documentation  of  other  validation  efforts.  A  broad  suite  of  vali¬ 
dation  tests  is  preferable  to  a  single  test  because  each  one  is  able  to  focus  on  a 
subset  of  processes  modeled  in  NUFT.  Therefore,  the  results  of  these  other  field 
validation  exercises  are  described  as  well. 


Validation  of  SVE  at  the  LLNL  B-518  Site 

A  detailed  description  of  the  SVE  validation  at  the  Lawrence  Livermore 
National  Laboratory  (LLNL)  Building  518  (B-518)  site  is  given  in  Nitao,  Martins, 
and  Ridley  (2000).  SVE  is  a  commonly  used  remediation  method  for  removing 
volatile  organic  contaminants  (VOCs)  from  the  vadose  zone.  Soil  gas  is  extracted 
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from  the  subsurface  by  appl5dng  a  vacuum  to  one  or  more  boreholes.  In  many 
cases,  clean  air  is  also  injected  at  other  locations  to  enhance  air  flow  rates. 

Building  518  was  constructed  in  1958  and  used  as  a  gas  cylinder,  solvent 
drum,  and  oil  drum  facility.  Several  sites  around  the  former  facility  were  identi¬ 
fied  to  have  potentially  high  concentrations  of  VOC,  primarily  TCE.  In  1993  a 
2-day  field  test  was  performed  at  a  borehole  located  in  the  area  of  the  highest  soil 
concentrations  to  demonstrate  the  treatability  of  the  site.  The  mean  permeability  of 
the  model  was  calibrated  based  on  the  wellhead  pressure  and  the  total  extracted 
vapor  flux.  The  initial  total  contaminant  mass  was  calibrated  using  the  vapor 
concentrations  from  the  borehole.  The  shape  of  the  vapor  concentration  history 
curve  agreed  well  with  calibrated  model  predictions  (Figure  13).  It  was  found  that 
the  initial  mass  had  to  be  increased  up  to  five  times  the  initial  estimates  using 
simple  spatial  interpolation,  indicating  that  initial  estimates  were  probably  too 
low. 


In  September  1995,  actual  remediation  of  the  site  was  begun  using  vapor 
extraction  from  the  same  borehole.  The  remediation  was  modeled  using  the  model 
calibrated  in  1993  with  the  focus  on  the  first  19  months  of  extraction  because  after 
that  period  other  extraction  boreholes  began  operation,  which  would  extend  the 
range  of  remediation  beyond  that  of  the  1993  calibration.  The  total  vapor  flux 
history  measured  in  the  field  was  input  to  the  calibrated  NUFT  model  in  the  form 
of  a  specified  flux  well  condition.  The  resulting  comparison  with  the  field  data 
and  model  prediction  is  shown  in  Figure  14.  Increasing  the  initial  mass  in  the 
model  by  10  percent  shows  the  sensitivity  to  the  mass  estimate.  Most  of  the  VOC 
stream  is  TCE,  but  other  relatively  minor  volatile  contaminants  that  are  present, 
but  not  included  in  the  model,  will  cause  total  VOC  to  be  underpredicted. 
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Figure  14.  Comparison  of  cumulative  total  VOC  mass  from  field  measurements  of  TCE  beginning 

September  1995  versus  those  predicted  by  the  model  during  first  19  months  of  SVE,  1993 

Using  the  same  homogeneous  model  that  was  calibrated  against  the  2-day 
treatability  test  in  1993,  very  good  agreement  with  the  produced  stream  was 
obtained  during  the  first  19  months  of  remediation  using  SVE  from  the  same  well. 
Although  the  main  goal  of  the  study  was  model  validation  against  SVE,  the  study 
also  showed  how  a  numerical  model  can  be  used  to  improve  initial  contaminant 
mass  estimates. 


Modeling  of  bioventing  at  Site  280,  Hill  AFB 

An  initial  modeling  effort  was  performed  for  bioventing  at  Site  280,  Hill 
AFB,  near  Ogden,  UT.  Bioventing  is  the  remediation  of  hydrocarbons  from  the 
vadose  zone  by  injecting  air  into  the  vadose  zone  to  enhance  aerobic  biodegrada¬ 
tion  of  hydrocarbons  from  naturally  occurring  microbes.  Compared  to  SVE,  the 
required  flow  rates  are  low.  The  vadose  zone  at  Site  280  was  contaminated  with 
jet  fuel. 

The  components  modeled  were  water,  N2,  O2,  CO2,  and  benzene  (CeHe).  All 
components  can  partition  into  both  gaseous  and  aqueous  phases,  which,  for  exam¬ 
ple,  allows  for  evaporation  of  water  and  dissolution  of  benzene  in  the  aqueous 
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phase.  Both  diffusive  and  advective  transport  of  all  components  in  both  fluid 
phases  were  considered.  Since  no  NAPL  has  been  observed  at  the  site,  no  free 
phase  was  considered. 

Stoichiometry  of  the  biodegradation  process  was  described  by  the  following 
reaction 


CfiHe  +  7.5  O2  =>  6  CO2  +  3  H2O  (1) 

An  irreversible,  dual-substrate  Monod  kinetic  reaction  was  used  to  model  this 
reaction.  The  coefficient  for  the  rate  equation  was  cahbrated  based  on  plume  size 
and  total  mass  removal. 

Modeling  without  bioventing  showed,  as  expected,  that  biodegradation  is 
limited  by  the  lack  of  oxygen.  Natural  biodegradation  was  predicted  to  destroy 
15.4  percent  of  the  contaminant  mass  while  bioventing  was  predicted  to  destroy 
99.6  percent.  Further  work  ended  because  of  re-prioritization  from  the  SERDP 
Office.  Sun  et  al.  (2000)  documents  this  study. 


Modeling  dynamic  steam  stripping— the  Visaiia  Test 

In  the  summer  of  1997  a  test  was  performed  at  Visalia,  CA,  that  injected 
steam  into  a  creosote-contaminated  formation  (Newmark  and  Aines  1998; 
Newmark  et  al.  1999).  Water  sampling  and  temperature  monitoring  were  per¬ 
formed  at  an  extraction  well  23  m  (75  ft)  from  the  extraction  well.  A  NUFT  model 
was  developed  that  incorporated  heterogeneous  distribution  of  hydrological  prop¬ 
erties  as  given  in  Figure  11.  The  predicted  breakthrough  of  the  steam  front  at  the 
extraction  well  was  approximately  48  hours,  which  was  very  close  to  the  observed 
breakthrough  time. 

Several  noble  gas  tracers  added  to  injected  water  and  steam  to  track  their 
movement  were  subsequently  monitored  at  the  extraction  well.  Figure  15  shows 
the  breakthrough  curves  at  the  extraction  well  from  xenon  that  was  injected  along 
with  the  steam.  The  different  curves  represent  different  initial  conditions.  The 
decrease  in  measured  concentrations,  as  shown  in  the  figure,  that  begins  at 
1 5  hours  was  caused  by  partial  failure  of  the  extraction  well  due  to  emulsification 
of  NAPL.  The  NUFT  model  did  not  include  the  pump  failure  because  its  exact 
degree  of  impact  is  uncertain.  Also,  it  should  be  noted  that  most  of  the  simulations 
were  done  before  the  field  test  was  actually  performed.  The  agreement  between 
field  data  and  model  prediction  is  considered  to  be  good  considering  the  complex¬ 
ities  involved  in  the  steam  stripping  process  and  the  high  degree  of  rmcertainties 
involved  in  modeling  a  site  that  had  relatively  little  quantitative  characterization  of 
hydrologic  properties. 


Benchmarking  of  NAPL  movement  and  steam  stripping 

NUFT  has  the  capability  to  model  heating  of  porous  media  through  AC  elec¬ 
trical  heating  (Carrigan  and  Nitao  2000).  LLNL  has  been  performing  scoping 
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Figure  15.  Breakthrough  curves  for  xenon  tracer  that  was  injected  with  the  steam 
(from  Newmark  and  Aines  1998) 

studies  to  model  the  enhancement  of  production  from  heavy  oil  reservoirs  by  this 
process  (Sahni,  Kxunar,  and  Knapp  2000).  The  sponsor  requested  that  NUFT  be 
verified  against  an  established  petroleum  industry  reservoir  simulator,  in 
particular,  the  CHEARS  code  developed  by  Chevron  Oil  Company  (Chien  and 
Northrup  1993).  The  selected  test  problem  was  the  steam  flooding  of  a  hypo¬ 
thetical  reservoir.  Of  course,  simulation  of  this  condition  has  implications  for 
modeling  recovery  of  NAPL  source  contamination. 

CHEARS  was  run  in  this  exercise  as  a  black  oil  simulator;  that  is,  the  petro¬ 
leum  is  treated  as  a  single  nonvolatile  component.  Because  the  NUFT-USNT 
module  cannot  model  a  nonvolatile  component,  the  oil  phase  had  small  amounts 
of  dissolved  noncondensable  gas.  The  benchmark  problem  consisted  of  a  three- 
dimensional  rectangular  domain  (109  m  (358  ft)  wide  by  109  m  (358  ft)  deep  and 
18.3  m  (60.04  ft)  tall)  with  a  steam  injection  well  at  the  center  and  a  fully  pene¬ 
trating  production  well  at  each  comer  of  the  domain.  The  formation  was  a  perme¬ 
able  sandstone  with  two  low-permeability  shale  zones  of  thickness  1.524  m  (5  ft), 
starting  at  3.05  m  (10  ft)  from  the  bottom  of  the  model  and  3.05  m  (10  ft)  from  the 
top  of  the  model.  Steam  was  injected  only  in  the  sandstone  between  the  two  shale 
zones.  The  bottom  and  top  of  the  model  formed  impervious  boundaries.  The 
formation  was  tilted  by  10  degrees  from  the  vertical  along  one  of  the  sides  of  the 
model. 

Figure  16  shows  the  oil  phase  (NAPL)  saturation  at  time  0.58  year  from  the 
start  of  the  simulation  for  a  cross  section  that  runs  through  the  injection  well.  The 
formation  is  tilted  upwards  going  from  the  left  to  right.  An  oil  bank,  as  indicated 
by  the  red,  developed  around  the  steam  injection  well.  The  right  portion  of  the 
bank  that  is  shown  progressed  farther  because  of  the  imposed  tilt.  The  bluish 
zones  are  the  shale  zones,  which  have  very  little  oil  saturation.  Figure  17  shows 
the  plan  view  of  the  oil  saturation  at  the  same  time  period,  at  a  horizontal  plane 
underneath  the  upper  shale  zone,  clearly  showing  the  resulting  oil  bank. 
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Figure  16.  Oil  phase  saturation  along  a  vertical  cross  section  through  the 

injection  well  parallel  to  the  sides  of  the  model  (the  model  tilts  upward 
by  1 0  degrees  from  left  to  right) 


Figure  17.  Plan  view  of  oil  phase  saturation  in  the  sand  layer  Just  below  the  upper 
shale  layer  (the  model  tilts  upwards  from  left  to  right) 


Figure  18  shows  the  total  oil  flux  and  water  flux  produced  from  the  field  as 
predicted  by  the  CHEARS  and  NUFT  models.  The  underprediction  of  oil  by 
NUFT  compared  with  the  CHEARS  model  is  probably  due  to  the  reduction  in  oil 
flow  from  evolution  of  noncondensable  gas  around  production  wells  in  the  NUFT 
calculation.  Comparison  between  the  models  is  considered  good,  especially  when 
compared  to  other  industry  code  benchmarking  exercises  (Aziz,  Ramesh,  and 
Woo  1987).  Note  that  differences  can  also  arise  from  the  manner  in  which  the 
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Figure  18.  Comparison  of  oil  and  water  fluxes  from  the  NUFT  and  CHEARS 
models  (courtesy  of  R.  Knapp,  Geosciences  and  Environment 
Technologies  Division,  Lawrence  Livermore  National  Laboratory, 
Livermore,  CA,  6  October  2000) 


models  control  their  time-steps  and  other  aspects  of  their  numerical  solution 
algorithms. 


Implementation  into  the  GMS  Interface 

A  sophisticated  NUFT  interface  for  GMS  was  developed  using  information 
about  NUFT  input  from  the  author  of  NUFT.  ‘  Because  most  environmental 
remediation  methods  modeled  by  NUFT  use  the  USNT  module,  it  was  decided  to 
implement  only  this  module.  Because  of  the  many  options  available  in  NUFT  and 
because  of  its  generality,  considerable  programming  effort  was  required  to 
program  the  GMS-NUFT  interface. 

Some  modifications  to  the  NUFT  code  were  also  required  because  NUFT  was 
originally  written  to  be  run  at  the  MS-DOS  or  Unix  command  line,  instead  of 
interactively  under  a  graphical  user  interface,  such  as  GMS.  First,  a  new  output 
option  in  NUFT  was  implemented  and  tested  that  allows  NUFT  to  create  output 
files  in  the  standard  GMS  binary  output  format  so  that  GMS  can  read  NUFT 
output  directly  for  graphical  output.  The  output  times  for  the  data  can  be  specified 
by  the  user  through  the  GMS  interface.  Next,  NUFT  and  its  program  build  system 
were  ported  from  the  GNU  C++  compiler  to  the  Microsoft  Visual  C++  compiler 


‘  Personal  Communication,  1999,  J.  J.  Nitao,  Lawrence  Livermore  National  Laboratory, 
Livermore,  CA. 
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to  support  long  file  names  under  Windows  NT.  NUFT  assumes  that  all  user 
output  files  should  be  created  in  the  directory  in  which  it  is  being  run  whereas 
GMS  passes  the  full  directory  path  of  the  input  file  and  expects  all  output  files  to 
go  in  the  directory  of  the  input  file.  A  function  supplied  by  Brigham  Young 
University  (BYU)  to  extract  the  executable  and  input  file  directory  paths  was 
integrated  into  NUFT  to  be  consistent  with  the  GMS  convention.  A  similar 
change  was  made  with  regard  to  the  directory  location  of  the  NUFT  start-up  files, 
which  is  now  based  on  parsing  the  full  directory  path  of  the  NUFT  executable.  All 
of  these  changes  had  to  be  made  while  maintaining  compatibility  with  non-PC 
platforms. 

The  following  steps  for  building  a  NUFT  simulation  using  GMS  are  from  the 
GMS  help  utility: 

a.  A  cell-centered  three-dimensional  grid  is  constructed. 

b.  The  basic  NUFT  options  are  selected  including  general  options, 
equations,  time-steps,  and  solver. 

c.  Regions  of  the  grid  are  selected  and  marked  as  ranges. 

d.  All  boundary  conditions  and  sources/sinks  are  assigned  to  ranges.  The 
ranges  also  define  material  zones. 

e.  The  material,  phase,  and  component  properties  are  defined. 

/  Boundary  conditions  and  phase  and  component  sources  are  assigned  to 
predefined  ranges. 

g.  Wells  are  created  and  assigned  the  appropriate  phase  and  component 
fluxes. 

h.  Initial  conditions  are  defined. 

i.  Output  control  options  are  selected. 

j.  The  model  is  saved  and  NUFT  is  launched  firom  the  GMS  menu. 

k.  The  NUFT  solution  is  read  into  GMS  for  postprocessing. 

An  existing  NUFT  simulation  can  be  read  into  GMS  using  the  “Read  Simulation” 
command  in  the  NUFT  menu.  Once  a  NUFT  simulation  is  saved  to  disk,  NUFT 
can  be  launched  by  selecting  the  “Run  NUFT”  command  firom  the  NUFT  menu. 
The  command  brings  up  a  dialog  listing  the  path  to  the  NUFT  executable  and  the 
most  recently  saved  NUFT  simulations.  In  most  cases,  these  paths  do  not  need  to 
be  edited.  When  the  “OK”  button  is  selected,  NUFT  is  launched  and  the  input  file 
is  passed  to  NUFT  as  a  command  line  argument.  NUFT  is  launched  in  a  separate 
window  and  the  console  output  fi'om  NUFT  is  displayed. 

Figure  19  shows  an  example  dialog  menu  for  specifying  the  phase-  dependent 
material  properties  of  a  NUFT  material  type. 
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Figure  19.  GMS  dialog  menu  in  the  NUFT  interface  for  specifying  the 

phase-dependent  material  properties  of  a  single  material  type 

Recommendations 

Because  multiple  flow  and  transport  mechanisms  occur  during  most  remedia¬ 
tion  processes,  preferably  more  than  one  field  test  should  be  used  to  verify  compu¬ 
ter  models.  Each  test  in  a  broad  suite  of  tests  can  be  used  to  focus  the  vaUdation 
on  a  particular  subset  of  dominant  processes. 

Table  4  is  a  matrix  diagram  listing  various  physical  processes  modeled  by 
NUFT  with  x’s  indicating  appropriate  tests  validating  each  process.  If  a  process 
plays  only  a  minor  role  in  a  test,  it  is  not  checked.  For  example,  aqueous  phase 
transport  due  to  infiltration  occurs  in  the  LLNL  B-518  SVE  remediation  test,  but 
at  that  site  it  is  unimportant  relative  to  gaseous  phase  transport,  so  it  is  not 
marked.  However,  in  less  arid  regions,  aqueous  phase  transport  can  play  a  sig¬ 
nificant  role  in  contaminant  transport  during  SVE,  at  least  during  periods  of  high 
raiirfall.  Aqueous  phase  transport  is  an  important  part  of  SVE  in  such  circum¬ 
stances.  Thus,  it  is  important  to  have  other  tests  for  validating  aqueous  phase 
advection  such  as  a  trench  infiltration  test,  which  has  been  performed. 

The  Visalia  field  test  and  the  steam-flood  simulator  test  are  important  for 
modeling  thermal  remediation  processes  and  for  problems  where  there  is  NAPE 
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Table  4 

NUFT  Process  Validation  Matrix 

Main  Processes  Validated 

Field  Tests  || 

A 

B 

c 

D 

E 

F 

Gaseous  flow 

X 

X 

X 

X 

Gaseous  transport 

X 

X 

X 

Aqueous  flow 

X 

X 

X 

X 

Aqueous  transport 

X 

X 

X 

X 

NAPL  flow 

X 

Steam  flow 

X 

X 

Steam  transport 

X 

Thermal  advective  transport 

X 

X 

X 

Thermal  conductive  transport 

X 

X 

X 

Fluid  phase  changes 

X 

X 

Heterogeneous  systems 

X 

X 

Barometric  pumping  effects 

X 

Transport  in  fractures 

X 

Infiltration  flow 

X 

Infiltration  transport 

X 

Injecting  and  producing  wells 

X 

1 

X 

X 

Notes: 

Field  Tests 

A.  LLNL  B-518  SVE  remediation 

B.  Nuclear  treaty  test  (Carrigan  et  al.  1996) 

C.  CAMBRIC  trench  study,  personal  communication  with  F.  B.  Thompson,  Geosciences  and 
Environmental  Technologies  Division,  Lawrence  Livermore  National  Laboratory,  15  November 
2000. 

D.  Cheshire  nuclear  test  (Maxwell  et  al.  2000;  Tompson  et  al.  2000) 

E.  Visalia  field  test 

F.  Steam-flood  simulator  benchmark 

X  =  Appropriate  test  for  validating  the  process. 

movement.  NAPL  movement  is  important  at  SVE  sites  where  there  is  free 
product. 

The  results  of  the  validation  tests  modeled  by  NUFT  were  very  encouraging. 
Although  it  is  impossible  to  rigorously  prove  the  validity  of  NUFT,  or  any  com¬ 
puter  code  for  that  matter,  the  tests  that  were  performed  give  confidence  that  the 
code  can  be  used  for  modeling  SVE,  thermal  remediation  methods  such  as 
dynamic  steam  stripping,  and  general  prediction  of  vadose  zone  flow  and  trans¬ 
port  processes.  However,  it  should  be  recognized  that  NUFT  is  a  fairly  complex 
simulator  offering  much  flexibility.  Thus,  it  will  take  time  for  this  simulator  to 
arrive  at  a  state  where  it  is  widely  applied. 
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Model  Description 

OS3D  is  a  simulator  for  multidimensional,  multicomponent  reactive  transport 
in  porous  media.  It  predicts  the  behavior  of  multiple  interacting  chemical  compo¬ 
nents  by  modeling  subsurface  reaction  and  transport  processes  in  a  quantitatively 
mechanistic  framework  that  accounts  for  spatially  variable  geochemical  conditions 
and  properties.  These  reactive  components  include  minerals,  aqueous  chemical 
species,  and  species  sorbed  to  mineral  surfaces.  The  modeling  technology  spans 
multiple  disciplines  (e.g.,  hydrology,  geology,  chemistry,  biology,  numerical 
analysis,  and  computational  methods)  to  address  multiple  end  user  applications 
(e.g.,  diagenesis,  mineral  extraction,  remediation,  and  waste  management). 

OS3D,  originally  developed  by  Steefel  and  Yabusaki  (1996)  in  1995,  has 
been  updated  over  the  intervening  years.  This  latest  version  of  OS3D  is  written 
entirely  in  the  FORTRAN  90  language  and  is  now  part  of  a  reactive  transport 
package  called  CRUNCH  (Steefel,  in  preparation).  The  GMS  implementation  of 
OS3D  and  the  development  of  a  parallel  processing  version  is  a  collaborative 
effort  among  ERDC,  LLNL,  and  Pacific  Northwest  National  Laboratory. 

Transient  and  steady-state  phenomena  can  be  addressed  in  batch,  reaction 
path,  one-,  two-,  and  three-dimensional  systems  that  include  advective,  dispersive, 
and  diffusive  transport.  The  OS3D  reactive  transport  simulator  can  be  used  in 
conjunction  with  nonuniform,  single-phase-fluid  velocity  fields  provided  by  the 
MODFLOW  (McDonald  and  Harbaugh  1988)  groundwater  flow  simulator.  OS3D 
solves  a  set  of  nonlinear  partial  differential  equations  describing  coupled  reactions 
and  transport  m  a  multicomponent  system  using  a  finite  volume  discretization  of 
the  problem  domain.  A  mixed  system  of  equilibrium  and  kinetic  reactions  repre¬ 
senting  aqueous  and  siuface  complexation,  ion  exchange,  precipitation,  and 
dissolution  can  be  solved  (Steefel  2000).  An  extended  Debye-Huckel  model  is 
used  to  calculate  activity  corrections.  OS3D  uses  an  extensible  version  of  the 
EQ3/EQ6  (Wolery  1992)  temperature-dependent  thermodynamic  database  for 
reaction  stoichiometry,  stability  constants,  mineral  molar  volume,  and  activity 
model  parameters.  Constant  temperature  and  temperature  gradient  conditions  can 
be  simulated. 

Total  component  concentrations  are  transported  in  both  the  dispersion  and 
advection  routines,  and  form  the  basis  for  mass  balance  equations  in  the  reactive 
chemistry  step.  The  actual  unknowns,  however,  are  the  individual  “primary”  or 
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“component”  species  concentrations.  This  approach  effectively  folds  a  speciation 
calculation  into  the  nonlinear  ordinary  differential  equations  describing  the  con¬ 
servation  of  solute  mass  in  the  system.  The  approach  also  allows  for  mixed 
equilibrium  and  kinetic  reactions  within  the  overall  reaction  network.  The 
following  are  major  reactive  chemistry  components  included  in  OS3D: 

a.  Equilibrium  aqueous  complexation. 

b.  Kinetically  controlled  aqueous  phase  reactions  (e.g.,  redox  and 
biological). 

c.  Mineral  dissolution/precipitation  kinetics. 

d  Multicomponent  ion  exchange. 

e.  Surface  complexation  with  a  double-layer  model  used  for  an  electrostatic 
correction. 

f.  Redox  catalysis  on  mineral  surfaces. 

g.  Radioactive  decay. 

h.  Linear  adsorption  isotherms  Kj. 

In  theory,  there  are  no  restrictions  on  the  number  of  different  processes 
included  in  a  simulation  since  physical  memory  for  each  of  these  is  allocated  at 
run  time.  The  code  is  also  accompanied  by  an  initialization  routine  that  allows  the 
use  of  a  variety  of  different  constraints  (e.g.,  total  concentration,  equilibrium  with 
a  mineral  phase,  equilibrium  with  a  gas  phase,  charge  balance)  to  speciate  initial 
or  boundary  conditions.  Given  a  choice  of  primary  component  species,  the  code 
searches  a  version  of  the  EQ3/EQ6  database,  which  has  been  augmented  with 
surface  complexation  and  ion  exchange  thermodynamic  data  along  with  kinetic 
rate  laws  for  aqueous  and  mineral  reactions. 

At  startup,  the  code  allows  the  specification  of  a  number  of  different  zones 
within  the  domain  of  interest,  which  have  distinct  initial  solute  concentrations, 
mineral  abimdances,  and/or  grid  spacing.  The  code  also  includes  a  distribution  of 
species  calculation  to  specify  the  initial  and  boundary  solute  concentrations.  The 
code  reads  a  thermodynamic  database  based  on  the  EQ3/EQ6  database,  which 
covers  temperatures  between  0  and  300  °C  along  the  water  saturation  curve  for  a 
large  number  of  mineral  and  aqueous  species.  The  code  also  allows  for  a  range  of 
mineral-water  reaction  rate  laws,  including  parallel  reactions  (e.g.,  pH-dependent 
and  pH-independent  rates)  and  proton-promoted  and/or  ligand-promoted 
dissolution. 

The  operator  splitting  solution  methodology  comprises  an  explicit,  third-order 
accurate  TVD  advection  scheme;  an  implicit,  second-order  accurate  dispersion 
and  diffusion  scheme;  and  an  implicit  backward  Euler  scheme  for  the  system  of 
algebraic  equations  representing  the  differential  equations  for  reactions.  Concen¬ 
tration  fronts  are  accurately  preserved  with  the  high-resolution,  shock-capturing 
TVD  method. 

An  adaptive  time-stepping  algorithm  is  used  to  control  the  convergence  of  the 
nonlinear  iteration  scheme  within  the  limits  of  the  maximum  allowable  time-step 
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that  maintains  the  stability  of  the  explicit  numerical  scheme  for  advective 
transport.  Fully  coupled  or  “global  implicit”  approaches  can  also  provide  solutions 
to  reactive  transport  problems;  however,  OS3D  has  some  real  advantages  in  terms 
of  memory  usage,  CPU  time  required  per  time-step,  the  ability  to  handle  up  to 
three-dimensional  transport,  and  the  coupling  to  higher  order  advective  transport 
schemes.  This  made  it  an  attractive  choice  for  incorporation  into  the  GMS 
software  suite. 


Inquiries  about  the  code  should  be  directed  to: 


Carl  I.  Steefel 

Lawrence  Livermore  National  Laboratory 
Geosciences  and  Environmental  Technologies 
Phone:  925-424-9807 
FAX:  925-423-1997 
E-mail:  steefel@llnl.gov 

Postal  Address:  Street  Address: 


LLNL,  L-204 
PO  Box  808 

Livermore,  CA  94551-9900 


LLNL,  L-204 
7000  East  Avenue 
Livermore,  CA  94550-9900 


Model  Enhancements 

OS3D  contains  a  number  of  enhancements  over  previous  versions  of  the  code 
that  have  been  released  (Steefel  and  Yabusaki  1996).  The  present  version  has 
been  entirely  rewritten  in  FORTRAN  90,  incorporating  data  stracture  modules, 
run  time  dimensioning  of  arrays,  and  dynamic  memory  allocation.  These  improve¬ 
ments  allow  (a)  problem  initialization  routines  to  be  executed  without  loading 
memory  for  the  global  spatial  arrays;  (b)  problem  restart  that  skips  the  initiali¬ 
zation  routines  altogether;  and  (c)  more  information  that  can  be  used  in  the 
parallel  implementation.  A  keyword-based  input  has  been  implemented  that 
allows  the  entry  of  data  categories  in  any  order  and  the  ability  to  comment  entries 
so  they  will  be  ignored  by  the  code. 

OS3D  now  allows  for  the  importation  of  a  transient  or  steady-state  flow  field 
from  MODFLOW.  With  this  added  capability  comes  the  ability  to  simulate  a 
variety  of  source  terms,  including  injection  and  pumping  wells,  drains,  rivers,  and 
other  fixed-head  boundary  conditions.  New  routines  were  developed  to  accom¬ 
modate  data  transfer  from  MODFLOW  and  the  GMS  interface  to  OS3D. 
MODFLOW  input  data  are  now  transferred  through  the  GMS  interface  into  the 
OS3D  input  file.  In  the  GMS,  MODFLOW  initial  and  boimdary  conditions  and 
flow  fields  are  written  to  output  files  for  use  by  transport  routines  in  the 
MT3DMS  (Zheng  and  Wang  1999)  simulator.  These  same  files  are  now  read  by 
OS3D.  This  condition  required  the  modification  of  the  OS3D  transport  routines  to 
be  consistent  with  the  MODFLOW  for  (a)  three-dimensional  cell  thickness 
variability;  (b)  flows  (rather  than  Dacy  velocities)  that  are  specified  at  cell  faces; 

(c)  source/sinks  (recharge,  evapotranspiration,  wells,  drains,  rivers,  streams);  and 

(d)  inactive  cells. 
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Enhancements  to  the  chemical  processes  simulated  with  OS3D  include  the 
following: 

a.  Option  to  automatically  include  every  relevant  secondary  chemical 
species  in  the  database  after  specifying  a  set  of  primary  species.  This  is 
typically  done  early  in  the  development  of  a  problem  specification  when  a 
geochemically  consistent  set  of  component  concentrations  is  being 
initialized  to  ensure  that  key  reactions  are  not  omitted. 

b.  Multicomponent  (multivalent)  ion  exchange  with  multiple  exchange  sites. 
The  user  may  choose  to  use  either  the  Gaines-Thomas,  Gapon,  or 
Vanselow  conventions  for  calculating  exchange  species  activities  (e.g., 
Appelo  andPostma  1993). 

c.  Multisite  surface  complexation,  with  or  without  an  electrostatic  correction 
accounting  for  the  effect  of  surface  charge  on  mineral  surfaces.  The 
electrostatic  or  “coulombic”  correction  is  based  on  the  double  layer  model 
(Davis  and  Kent  1990;  Parkhurst  1995;  Bethke  1996). 

d.  Aqueous  phase  (homogeneous)  kinetics  including 

(1)  Biodegradation  using  Monod-type  rate  laws  or  variable  order 
(zeroth-,  first-,  second-order)  rate  laws. 

(2)  Radioactive  decay  chains. 

e.  Solid  phase  reactions  using  either  transition  state  theory  types  of  rate  laws 
with  various  catalytic  or  inhibitory  fimctions  or  with  Monod-type  kinetics. 

Unlike  previous  versions  of  the  code,  only  total  component  concentrations  are 
transported  in  the  present  version,  so  no  basis  switching  (Steefel  and  Yabusaki 
1996)  is  allowed. 


Model  Validation 

The  OS3D  simulator  was  applied  to  a  field  study  (Battelle  1998)  of  a  pilot- 
scale  fiinnel  and  gate  treatment  system  at  Moffett  Field,  California,  that  used  zero- 
valent  iron  (Fe°)  to  engineer  TCE  destruction  (Yabusaki  et  al.  2001).  To  date,  Fe® 
has  been  emplaced  in  a  variety  of  subsurface  barrier  geometries  in  a  variety  of 
subsurface  settings,  e.g.,  Moffett  Field,  Elizabeth  City,  Dover  AFB,  and  Denver 
Federal  Center  (U.S.  Environmental  Protection  Agency  1999).  While  it  is  clear 
that  Fe”  can  be  effective  in  transforming  chlorinated  hydrocarbons  to  benign 
byproducts,  there  are  many  other  components  in  the  groundwater  that  will  also  be 
affected  by  the  presence  of  Fe“.  Of  concern  are  reactions  that  could  adversely 
affect  the  long-term  performance  of  the  iron  barrier  or  create  secondary  contami¬ 
nants.  Since  the  solution  chemistry  (e.g.,  pH,  Eh,  alkalinity,  dissolved  gases,  major 
cations  and  anions)  of  groundwater  will  vary  fi:om  site  to  site,  the  long-term  per¬ 
formance  of  a  “standard”  implementation  of  an  iron  barrier  may  also  be  expected 
to  vary.  Understanding  the  impact  of  the  local  groundwater  chemistry  on  the  per¬ 
formance  of  the  iron  barrier  provides  opportunities  to  assess  the  feasibility  of  the 
technology  and  adapt  the  engineering  design  to  site-specific  conditions. 
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Problem  description 

The  funnel  and  gate  methodology  emplaces  impermeable  walls  in  a  saturated 
aquifer  to  “funnel”  contaminated  groundwater  through  one  or  more  “gates”  where 
passive  treatment  takes  place  (Starr  and  Cherry  1994).  Shallow  (<  15  m  (50  ft) 
below  grormd  surface)  aquifers  imderlain  by  a  continuous  bottom-confining  layer 
are  the  best  candidates  for  the  technology.  Sheet  piles  with  sealable  joints  are 
typically  used  to  construct  the  impermeable  walls,  whereas  the  gates  contain 
reactive  cells  filled  with  varying  amormts  and  geometries  of  granular  Fe®.  As 
contaminated  groundwater  passes  through  the  reactive  cell,  transformation  of  the 
target  contaminants  occurs. 

In  the  Moffett  Field  case,  the  gate  is  5.5  m  (18  ft)  deep,  3.0  m  (10  ft)  long, 
and  3.0  m  (10  ft)  wide  (Figure  20).  The  reactive  iron  cell  is  as  wide  and  deep  as 
the  gate  but  is  1.8  m  (6  ft)  long  with  0.6-m-  (2-ft-)  long  pea  gravel  zones  on  either 
end.  The  pea  gravel  has  a  porosity  of  0.33.  The  granular  iron  in  the  reactive  cell 
was  manufactured  by  Peerless  Metal  Powders,  Inc.,  in  the  -8  to  +40  mesh  particle 
size  range.  The  iron  was  characterized  with  a  packed  porosity  of  0.66  and  a  spe¬ 
cific  surface  area  of  1 .5  m^/gm,  which  are  typical  for  these  materials  (Johnson, 
Scherder,  and  Tratnyek  1996). 

After  installation  of  the  funnel  and  gate  in  April  1996,  a  very  comprehensive 
field  data  collection  effort  comprising  five  quarterly  monitoring  events  after  the 
installation  of  the  permeable  reactive  barrier  was  performed  through  October  1997 
(Battelle  1998).  Chlorinated  volatile  organic  compounds  (CVOCs^  hydrocarbons, 
major  ions,  dissolved  gases,  temperature.  Eh,  pH,  and  alkalinity  were  measured. 
Unless  otherwise  noted,  the  measured  data  reported  are  from  Battelle  (1998), 
which  also  includes  a  description  of  the  sampling  and  analytical  procedures. 


Model  specification 

Selection  of  aqueous  components.  Based  on  the  observed  changes  to  the 
groundwater  chemistry  from  the  aquifer  through  the  iron  zone,  the  following 
aqueous  species  were  modeled: 

a.  Organic  chemistry:  TCE,  c/5-1, 2-DCE,  CH4,  C2H4,  CF 

b.  Carbonate  chemistry:  H*^,  OH~,  HCO3",  C02(aq),  CO{^ 

c.  Major  metal  ions:  Ca^^,  Mg^^ 

d.  Redox  couples:  Fe^^/Fe^^  O2  (aq)/H20,  H20/H2(aq),  SO^'/HS' , 
NO’/NH; 

e.  Metal  complexes: 

FeOH^^ ,  Fe(OH)2 ,  Fe(OH)3  (aq),  Fe(OH):, ,  Fe(OH)^ ,  Fe(OH)2  (aq), 
Fe(OH)3 ,  FeC03  (aq),  FeHCO; ,  CaHCO; ,  CaCOj  (aq),  MgECO^ , 
MgCOj  (aq),  CaSO^  (aq),  MgSO^  (aq),  FeSO^  (aq) 
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As  a  first  approximation,  the  homogeneous  reactions  were  modeled  as  equilibrium 
reactions  using  the  stoichiometry  and  thermodynamics  from  the  EQ3/EQ6 
database  (Wolery  1992).  A  temperature  of  25  °C  was  assumed  for  the  reaction 
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thermodynamics,  which  was  appropriate  given  observed  temperatures  in  the 
20-22  °C  range. 

Selection  of  minerals.  Other  than  the  chemical  reduction  of  aqueous  com¬ 
ponents,  the  precipitation  of  secondary  minerals  was  assumed  to  be  the  primary 
mechanism  for  removing  calcium,  magnesium,  carbonate,  hydroxide,  and  sulfate 
from  solution  in  the  iron  cell.  Based  on  observed  concentrations  in  the  field  and 
the  analyses  of  the  iron  cores,  a  set  of  candidate  minerals  was  individually  tested 
for  geochemical  consistency  (i.e.,  under  or  at  saturation  in  the  background  aquifer 
and  over  saturation  in  the  iron  cell)  using  the  batch  equilibrium  mode  of  the 
OS3D  simulator.  For  example,  CaCOaCs)  is  indicated  as  a  potential  secondary 
mineral,  but  it  is  aragonite,  rather  than  calcite,  that  is  thermodynamically  con¬ 
sistent  with  the  observed  concentrations. 

Based  on  these  equilibrium  modeling  analyses,  the  following  secondary 
minerals  were  included  in  the  model:  ferrihydrite  (Fe(OH)3(am)),  Fe(OH)2(am), 
FeS(am),  brucite  [Mg(OH)2],  siderite  (FeCOa),  aragonite  (CaCOs),  and  green  rust 
[Fe4“Fe2’“(0H)i2][S043H20].  Both  macldnawite  (Fei+xS)  and  “amorphous  FeS” 
have  nearly  the  same  free  energy  of  formation  (Lennie  and  Vaughan  1996). 
Although  subsequent  sulfidation  to  form  marcasite  and  pyrite  is  expected,  the 
exact  sulfidation  sequence  and  kinetics  of  each  step  have  not  been  characterized 
well  enough  to  be  included  in  the  model.  The  green  rust  reaction  thermodynamics 
are  from  Hansen,  Borggaard,  and  Sorensen  (1994).  The  reaction  thermodynamics 
for  the  remaining  secondary  minerals  are  from  the  25  °C  EQ3/EQ6  database. 

Reaction  rates.  The  degradation  kinetics  of  TCE  and  c/5-1, 2-DCE  were 
modeled  using  the  beta  elimination  pathway  (Roberts  et  al.  1996)  as  a  bulk  first- 
order  kinetic  reaction  that  includes  the  oxidation  of  iron  in  addition  to  the  trans¬ 
formation  of  the  chlorinated  organics  to  undifferentiated  hydrocarbons  and 
chloride.  In  this  case,  ethene  is  used  as  the  nominal  hydrocarbon  byproduct  as 
follows: 


3Fe“  +  TCE  +  3H^=>3Fe^^  +  3Cr  +  C2H4  (2) 

2Fe®  + c/5-1, 2-DCE  +  2H^  =>  2Fe+^  +  2Cr  +  C2H4  (3) 

Bench  scale  column  studies  (PRC  Environmental  Management,  Inc.,  1995)  using 
granular  iron  from  Peerless  Metal  Powders,  Inc.,  were  performed  on  50:50  mix¬ 
tures  of  iron  and  sand  to  determine  degradation  rates  for  these  two  reactions.  The 
results  were  scaled  to  the  100  percent  iron  and  packing  conditions  used  in  the 
field  resulting  in  TCE  half-life  of  0.4  hour  and  c/5-l  ,2-DCE  half-life  of 
1.35  hours. 

Other  than  the  iron  dissolution  that  occurs  in  conjunction  with  the  transfor¬ 
mation  of  the  chlorinated  hydrocarbons,  a  reversible  rate  law  for  mineral  precipi¬ 
tation  and  dissolution  is  used  that  is  based  loosely  on  transition  state  theory  (e.g., 
Lasaga  1984): 
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(4) 


-1 


where 

=  rate  of  precipitation  (>0)  or  dissolution  (<0)  of  mineral  m 
A„  =  surface  area  of  the  reacting  mineral 
k„  =  rate  constant 
Ns  =  number  of  species 
Oi  =  activity  of  a  species 
p  =  empirically  determined  exponent 
Q„  =  ion  activity  product 

K„  =  equilibrium  constant  for  the  mineral  -  water  reaction 

For  dissolution,  the  intrinsic  rate  constant  k„  is  multiplied  by  the  bulk  surface  area 
A„  (i.e.,  4.0  X  10^  mW  for  iron).  Precipitation  rates  in  this  study  are  reported  in 
units  of  mol/mVs  (Am  km). 

In  addition  to  the  iron  oxidations  in  Equations  2  and  3,  there  are  two  other 
parallel  iron  reactions  involving  the  reduction  of  dissolved  oxygen  and  reduction 
of  water: 


2Fe®  +  2H2O  +  O2  (aq)  <=>  2Fe^^  +  40H‘  (5) 

Fe®  +  2H2O  o  Fe^^  +  H2(aq)  +  20H'  (6) 

The  rate  law  (Equation  4)  for  the  dissolution  reaction  with  oxygen  (Equation  5) 
has  a  relatively  large  rate  constant,  4  x  10'^  moFm^/s,  which  is  multiplied  by  the 
02(aq)  activity.  It  is  the  only  inorganic  rate  law  in  the  model  that  has  a  first-order 
dependency  on  the  activity  of  a  species  (i.e.,  02(aq)).  The  rate  constant  for  the 
parallel  iron  dissolution  reaction  involving  the  reduction  of  water  (Equation  6) 
was  derived  fi-om  a  column  study  (PRC  Environmental  Management,  Inc.,  1995) 
using  the  granular  iron  from  Peerless  Metal  Powders,  Inc.,  and  Moffett  Field 
groundwater.  In  this  study,  a  I.22-m-  (4-ft-)  long,  10.16-cm-  (4-in.-)  diameter 
column  was  packed  to  a  porosity  of  0.38  with  a  50:50  mass  ratio  of  sand  to  iron. 
The  resulting  bulk  iron  surface  area  was  1.88  x  10®  m^W.  For  the  5-ml/min  feed 
rate  of  Moffett  Field  groundwater  (0.888  m/day  Darcy  flux),  rate  constants  for 
hydrolysis  ranging  fi-om  3  x  10  ‘^  to  5  x  10'*^  moFm^/s  best  captured  the  observed 
chemistry. 
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Table  5 

Mineral  Reaction  Rates 

Mineral 

Logk 

moles/m^/s 

Fe{OH)2(s) 

-5.9 

Ferrihydrite 

-5.9 

Siderite 

-6.9 

Aragonite 

-6.1 

FeS(am) 

-6.9 

Brucite 

-7.4 

Green  Rust 

-6.9 

Rate  laws  for  the  seven  secondary 
minerals  are  based  on  reversible 
mineral  reaction  rates  (Table  5)  from 
an  earlier  study  (Mayer  1999)  of  a 
reactive  iron  permeable  barrier  at  the 
U.S.  Coast  Guard  Support  Center  near 
Elizabeth  City,  NC.  Their  respective 
rates  were  adopted  for  the  calcium, 
magnesium,  and  sulfur  sinks  to  the 
aragonite,  brucite,  and  FeS/green  mst 
minerals  used  in  this  model. 


Model  geometry.  The  model  uses  a  one-dimensional  representation  of  the 
gate  that  includes  the  up-gradient  pea  gravel  zone  (0.6  m  (2  ft))  and  the  reactive 
iron  cell  (1.8  m  (6  ft)).  The  down-gradient  pea  gravel  zone  was  neglected  because 
of  contamination  from  the  down-gradient  aquifer.  The  system  is  discretized  into 
60  grid  cells  of  uniform  size  (0.04  m  (0.13  ft)). 

Flow  and  transport.  Pore  velocity  in  the  iron  cell  was  measured  by  a  variety 
of  techniques:  water  levels  (0.46-0.52  m/day  (1.5-1. 7  ft/day)),  tracer  tests 
(0.06-0.18  m/day  (0.19-0.59  ft/day)),  and  down-hole  velocity  meter  (0.34— 

1.86  m/day  (1.12-6.1  ft/day))  The  estimated  range  of  likely  pore  velocities  was 
0.06-0.6  m/day  (0.19-1.97  ft/day),  which  is  equivalent  to  a  Darcy  velocity  range 
of  0.04— 0.4  m/day  (0.13-1.3  ft/day)  (based  on  iron  porosity  of  0.66).  Dispersion 
and  diffusion  were  assumed  to  be  negligible. 

Boundary  and  initial  conditions.  The  boundary  condition  and  the  initial 
conditions  for  the  model  are  based  on  an  April  1997  sampling  event.  The  field 
data  are  used  to  create  two  sets  of  equilibrated  species  concentrations.  The  first  set 
is  used  for  the  up-gradient  boundary  condition  and  the  initial  condition  for  the  pea 
gravel,  whereas  the  second  set  is  used  for  the  initial  condition  of  the  iron  zone. 

The  significance  of  equilibrating  the  two  geochemical  conditions  is  based  on  the 
assumption  that  the  observed  conditions  are  at  quasi-equilibrium.  Thus,  there  was 
an  attempt  to  prevent  unrealistic  transient  behavior,  such  as  the  sudden  precipita¬ 
tion  or  dissolution  of  minerals  due  simply  to  the  initial  chemical  specification. 

A  geochemically  consistent  boundary  condition  was  generated  using  the 
OS3D  reactive  transport  simulator  that  reflects  the  measured  background  aquifer 
water  chemistry.  The  following  constraints  were  used: 

a.  pH  fixed  at  7. 1 ,  C02(aq)  at  4.0E-4. 

b.  Total  iron,  oxygen,  and  ethene  at  1  .OE-10  M,  reflecting  their  below- 
detection  concentrations. 

c.  Total  calciiun  at  3.9E-3  M. 

d.  Total  magnesium  at  2.7E-3  M. 

e.  Total  sulfate  at  3.5E-3  M. 

f.  Total  nitrate  at  3.9E-5  M. 
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g.  Total  chloride  at  1 .  lE-3  M. 

h.  TCEat9.1E-6M. 

/.  c/s-l,2,-DCE  at  2.2E-6  M. 


The  boundary  condition  was  also  the  initial  condition  for  the  pea  gravel  zone.  The 
porosity  of  this  zone  was  0.33  although  no  reactive  minerals  were  specified.  The 
boimdary  condition  concentrations  for  the  primary  components  are  summarized  in 
Table  6. 


1  Table  6 

1  Initial  and  Boundary  Conditions 

Species 

Boundary  Condition 
and  Pea  Gravel  Initial 

Condition 

(Molality) 

Initial  Condition 
for  Iron  Cell 
(Molality) 

7.94E-08 

3.16E-11 

5.93E-20 

1.16E-31 

3.07E-03 

7.50E-05 

Mg"" 

1.79E-03 

7.00E-05 

1  HCOa' 

2.28E-03 

3.55E-11 

Cl 

1.10E-03 

1.10E-03 

02(aq) 

1.00E-10 

2.26E-73 

TCE(aq) 

9.10E-06 

1.00E-10 

c/s-1 ,2-DCE(aq) 

2.20E-06 

4.20E-08 

C2H4(aq) 

1.00E-10 

7.80E-07 

d 

CO 

1.88E-03 

9.73E-11 

I  NOs- 

3.90E-05 

3.59E-82 

Mineralogically,  the  reactive  cell  is  initially  Fe®  occupying  34  percent  of  the 
volume.  A  geochemically  consistent  aqueous  initial  condition  for  the  reactive  cell 
was  calculated  by  OS3D  using  the  following  constraints:  pH  of  10.5;  an  arbitrarily 
small  concentration  (l.E-10  M)  specified  for  total  iron,  total  sulfate,  total  nitrate, 
total  carbonate,  total  oxygen,  and  TCE;  total  calcium  at  7.5E-5  M;  total  magnes¬ 
ium  at  7.0E-5;  total  chloride  at  l.lE-3  M;  ethene  at  7.8E-7M;  and  c/5'-l,2-DCE  at 
4.2E-8  M.  The  aqueous  initial  condition,  summarized  in  Table  6,  will  essentially 
be  flushed  out  of  the  system  after  1  pore  volume.  From  a  modeling  perspective, 
the  initial  condition  should  not  result  in  nonphysical  transients  that  significantly 
alter  immobile  phases. 


Results 

With  OaCaq)  nearly  absent  and  chlorinated  hydrocarbons  at  relatively  low  con¬ 
centrations  in  the  background  aquifer,  the  reduction  of  water  by  iron  oxidation  is 
the  dominant  iron  corrosion  reaction  as  well  as  the  principal  mechanism  for  ele¬ 
vating  pH  and  lowering  the  redox  potential  in  the  iron  cell  at  Moffett  Field.  pH 
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and  redox  conditions  contiol  the  behavior  of  many  of  the  chemical  components  in 
the  iron  cell;  consequently,  any  attempt  to  reproduce  spatially  variable  concen¬ 
trations  of  chemical  components  in  the  iron  cell  must  first  consider  the  interplay 
between  this  relatively  slow  hydrolysis  reaction  and  the  rate  of  transport,  which 
controls  the  residence  time  in  the  iron  cell.  Based  on  the  laboratory  column 
studies  and  the  field  observations,  there  is  a  range  of  possible  rates  for  the 
hydrolysis  reaction  (intrinsic  rate  constant  <  5  x  10'*^  mol/m^/s)  and  water  flux 
(0.04-0.40  m/day  (0.13-1.3  fit/day))  in  the  iron  cell.  In  this  case,  however,  there  is 
only  one  combination  that  can  systematically  accoimt  for  the  observed  behaviors 
in  the  iron  cell:  a  hydrolysis  reaction  rate  constant  of  5  x  10"*^  mol/m^/s  and  a 
Darcy  velocity  of  0.04  m/day.  A  slower  hydrolysis  reaction  rate  would  require  a 
Darcy  velocity  below  the  lowest  field  estimate,  whereas  a  faster  Darcy  velocity 
would  require  a  hydrolysis  reaction  rate  greater  than  observed  in  the  laboratory. 

The  subsequent  discussion  is  based  on  the  incorporation  of  the  aforemen¬ 
tioned  rates  into  an  OS3D  simulation  of  the  Moffett  Field  fimnel  and  gate 
treatment  system  for  one  year  of  operation,  April  1996  to  April  1997.  One¬ 
dimensional  modeling  results  are  compared  with  field  observations  from  14  mon¬ 
itoring  wells  in  the  gate,  of  which  8  are  designed  for  multilevel  sampling.  In 
general,  the  principal  variability  of  the  observed  aqueous  components  is  along  the 
longitudinal  axis  of  the  gate.  In  the  accompanying  figures,  the  modeling  results 
are  compared  with  the  projection  of  the  three-dimensional  field  data  onto  this  axis. 

Chlorinated  hydrocarbons.  The  modeled  TCE  concentrations  reflect  the 
rapid  degradation  observed  in  the  field  (Figure  21),  supporting,  to  some  degree, 
the  TCE  degradation  rate  derived  in  the  laboratory.  For  c/5-l,2-DCE,  the  labora¬ 
tory  rate  used  in  the  model  appears  to  be  high  (Figure  21)  as  the  cz5-l,2-DCE  is 
predicted  to  be  eliminated  after  0. 1  m  in  the  iron  cell  compared  with  0.5  m 
observed  in  the  field.  Relative  to  the  overall  system  of  chemical  components,  the 
chlorinated  hydrocarbon  chemistry  has  minimal  impact. 

pH.  The  modeled  pH  tracks  the  observed  rapid  increase  and  plateau  in  the 
iron  cell  very  well  (Figure  22).  The  hydrolysis  reaction,  which  acts  to  elevate  pH, 
is  always  active  throughout  the  iron  cell.  Conversely,  protons  are  continually 
released  during  the  precipitation  of  secondary  minerals,  buffering  the  initial  pH 
increase  at  the  iron  cell  entrance,  and  then  stabilizing  pH  at  10.4  for  the  final 
1.5  m  (5  ft). 

Eh-  Modeled  concentrations  of  redox  couples  were  used  to  calculate  Eh  under 
the  assumption  of  global  redox  equilibrium.  All  the  redox  couples  are  connected 
through  the  equilibrium  reactions  so  it  is  not  necessary  to  identify  the  dominant 
redox  couple.  The  calculation  is  based  on  half  redox  reactions; 


Eh 


23RT(  0 


(7) 
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concentration  (M) 


TCE  and  cis-1 ,2-DCE  versus  Distance:  1  Year 
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Figure  21.  Comparison  of  model  predictions  and  field  observations  after  1  year  of 
reactive  barrier  operation:  TCE  and  c/s-1 ,2-DCE 


pH  versus  Distance:  1  Year 


Distance  (m) 


Figure  22.  Comparison  of  model  predictions  and  field  observations  after  1  year  of 
reactive  barrier  operation:  pH 
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where 


R  =  universal  gas  constant 
T  =  temperature,  °K 

F  =  Faraday  constant,  electric  charge  of  1  mole  of  electrons,  96,500 
coulombs 

=  log  {equilibrium  constant) 


Ox  +  e'  =  Red 

Ox  =  oxygen  concentration 


For  25  °C,  2.3RT/F  is  0.059  v. 


The  model  reproduces  the  trend  of  positive  Eh  in  the  up-gradient  aquifer  and 
negative  Eh  in  the  iron  cell;  however,  all  predictions  are  more  extreme,  nearly 
double  the  field  measurements  (Figure  23).  The  reasons  for  this  lack  of  agreement 
cannot  be  determined  with  certainty.  However,  it  should  be  noted  that  measured 
species  concentrations  were  available  to  calculate  Eh  in  the  aquifer  for  two  of  the 
redox  couples.  In  this  case.  Eh  for  the  02(aq)/H20  (720  mV)  and  N2(aq)/N03' 

(690  mV)  redox  couples  compare  more  favorably  with  the  modeled  Eh  of  700  mV 
than  the  down-hole  probe  readings  of  ~300  mV. 


Eh  versus  Distance:  1  Year 


Distance  (m) 


Figure  23.  Comparison  of  model  predictions  and  field  observations  after  1  year  of 
reactive  barrier  operation:  Eh 
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Alkalinity  (as  CaCOs).  The  modeled  alkalinity  reproduces  the  steady 
decrease  from  0.003  M  to  0.0002  M  observed  in  the  first  0.5  m  (1.6  ft)  of  the  iron 
cell  (Figure  24).  The  decrease  in  alkalinity  in  the  iron  cell  is  due  primarily  to  the 
precipitation  of  carbonate  (siderite,  aragonite)  and  hydroxide  (Fe(OH)2,  ferrihy- 
drite,  brucite)  minerals.  In  the  model,  carbonates  are  completely  depleted  from 
solution  in  the  first  0.24  m  (0.79  ft)  of  the  iron  cell.  The  predicted  reduction  of 
aqueous  inorganic  carbonates  to  methane  is  a  relatively  minor  process,  having 
little  impact  on  alkalinity.  Furthermore,  the  highest  predicted  methane  concentra¬ 
tion  can  be  over  an  order  of  magnitude  smaller  than  methane  concentrations 
measured  in  the  reactive  cell.  The  elevated  methane  concentrations  indicate  the 
presence  of  an  uncharacterized  carbon  source  in  the  field;  carbon  impurities  in  the 
Peerless  iron  are  known  to  be  3-4  percent  (Fort  2000;  Phillips  et  al.  2000). 


Alkalinity  (as  CaC03)  versus  Distance:  1  Year 


Figure  24.  Comparison  of  model  predictions  and  field  observations  after  1  year  of 
reactive  barrier  operation:  alkalinity  (as  CaCOa) 


Dissolved  hydrogen.  H2(aq)  was  not  detected  (<4  x  10'^  M)  in  the  field 
despite  being  a  product  of  the  principal  iron  oxidation  reaction  that  drives  this 
chemical  system.  In  the  reactive  transport  model,  H2(ac[)  concentrations  increase 
tmtil  the  dissolved  gas  pressure  exceeds  the  hydrostatic  pressure  at  middepth  in 
the  iron  cell,  1.25  bars.  At  that  point,  hydrogen  gas  is  assumed  to  outgas,  main¬ 
taining  a  saturated  concentration  of  1  x  10'^  M. 

Major  ions.  In  the  reactive  transport  model,  magnesium  is  removed  from 
solution  by  the  precipitation  of  brucite.  The  predicted  magnesium  concentrations 
track  the  field  observations  very  well  (Figure  25),  matching  the  decline  and 
subsequent  leveling  from  2.7  x  10'^  M  to  less  than  1  x  10"*  M  over  the  first  0.3  m 
(10.98  ft)  of  the  iron  cell. 
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Figure  25.  Comparison  of  model  predictions  and  field  observations  after  1  year  of 
reactive  barrier  operation:  magnesium 

Sulfate  appears  to  have  similar  behavior  with  considerably  more  scatter  in  the 
field  measurements  (Figure  26).  Two  competing  pathways  for  the  removal  of 
sulfate  fi'om  solution  are  modeled:  precipitation  of  green  rust  and  reduction  of 
sulfate  to  sulfide  followed  by  precipitation  of  FeS(am).  While  the  precipitation  of 
green  rast  is  a  likely  sulfur  sink  in  the  iron  cell,  the  thermodynamics  of  this  reac¬ 
tion  limits  the  amount  of  sulfur  that  can  be  removed  from  solution.  In  contrast,  the 
precipitation  of  FeS(am)  is  capable  of  removing  essentially  all  of  the  aqueous 
sulfide  in  the  iron  cell.  Since  sulfate  concentrations  are  very  low  after  the  first 
0.5  m  (1.6  ft)  of  the  iron  cell  and  sulfide  is  always  below  detection,  it  appears  that 
aqueous  sulfate  that  could  not  be  precipitated  as  green  mst  is  rapidly  reduced  to 
sulfide  and  removed  from  solution  during  the  precipitation  of  FeS(am).  This  is 
precisely  how  sulfate  removal  in  this  system  is  modeled. 

In  the  iron  cell,  nitrate  transformation  is  so  rapid  and  complete  that  nitrate  is 
always  below  detection.  Nitrate  in  the  model  is  quickly  reduced  to  ammonia  via 
equilibrium  homogeneous  reactions  that  match  the  field  observations  quite  well 
(Figure  27).  N2(aq)  does  not  appear  to  be  part  of  the  nitrate/ammonia  redox 
chemistry  as  evidenced  by  fairly  unchanging  measurements  in  the  field. 

Calcium  is  also  rapidly  and  nearly  completely  eliminated  from  solution  in  the 
iron  zone.  In  this  case,  however,  the  precipitation  thermodynamics  of  the  sec¬ 
ondary  mineral,  aragonite,  can  accoimt  for  only  half  of  the  calcium  (~0.002  M) 
that  needs  to  be  removed  from  solution.  This  implies  that  additional  mechanisms 
for  calcium  to  be  removed  from  solution  need  to  be  investigated  (e.g.,  precipita¬ 
tion  of  other  calcium  minerals,  complexation  of  calcium  on  mineral  surfaces). 
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ure  26.  Comparison  of  model  predictions  and  field  observations  after  1  year  of 
reactive  barrier  operation:  sulfate 
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Figure  27.  Comparison  of  model  predictions  and  field  observations  after  1  year  of 
reactive  barrier  operation:  nitrate 
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Iron  dissolution.  The  hydrolysis  reaction  is  the  dominant  iron  dissolution 
mechanism.  In  the  first  0.3  m  (0.98  ft)  of  the  reactive  cell,  the  predicted  iron 
mineral  volume  fraction  is  reduced  by  0.44  percent.  Beyond  this  initial  zone,  the 
reduction  is  predicted  to  be  uniformly  0.38  percent. 

Secondary  minerals.  After  1  year  of  simulation,  the  precipitation  of  second¬ 
ary  minerals  is  predicted  to  be  highest  near  the  reactive  cell  entrance  where  arago¬ 
nite  and  siderite  dominate  the  new  mineral  assemblage  (Figure  28).  Porosity  is 
reduced  at  this  location  from  0.6600  to  0.6344.  FeS(am),  green  rust,  brucite,  and 
to  a  lesser  degree,  ferrihydrite,  are  actively  precipitating  over  the  initial  0.3  m 
(0.98  ft)  of  the  iron  cell,  with  moderate  accumulations  (mineral  volume  fractions 
<0.5  percent).  After  the  first  15  cm  of  iron,  the  carbonate  minerals  are  essentially 
absent  and  Fe(OH)2(am)  comes  to  dominate  the  rest  of  the  iron  zone  with  a 
1.4  percent  mineral  volume  fraction.  The  model  assumes  that  Fe(OH)2(am)  is  the 
stable  phase. 


Figure  28.  Predicted  mineral  volume  fractions  after  1  year  of  reactive  barrier 
operation 

Discussion 

Reactive  cell  performance.  The  Darcy  velocity  of  0.04  m/day  in  conjunction 
with  the  hydrolysis  rate  from  the  laboratory  (5  x  lO  '^  mol/mVs)  was  clearly  con¬ 
sistent  with  the  observed  behavior  of  pH,  alkalinity.  Mg,  SO^' ,  and  NO3 .  How¬ 
ever,  this  should  be  a  concern  because  the  preconstruction  Darcy  velocity  in  the 
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Mineral  precipitation  in  the  iron  cell  could  potentially  shorten  the  longevity  of 
the  reactive  barrier  and  thus  is  of  concern  to  the  future  deployment  of  similar  in 
situ  chemical  treatment  technologies.  In  the  initial  part  of  the  iron  zone  there  is  a 
net  reduction  in  porosity  of  0.015-0.03  each  year,  due  principally  to  the  precipita¬ 
tion  of  aragonite  and  siderite.  If  the  precipitation  is  confined  only  to  the  entrance 
of  the  iron  zone,  there  could  be  significant  reduction  in  porosity  and  permeability 
over  many  years  that  would  result  in  the  buildup  of  head  on  the  up-gradient  side 
of  the  iron  barrier.  Furthermore,  the  predicted  volume  of  out-gassed  hydrogen  is 
significant,  and  it  is  possible  that  hyi-ogen  gas  bubbles  are  occupying  a  portion  of 
the  void  space,  resulting  in  a  reduction  of  permeability.  Precipitation  of  secondary 
minerals  also  has  the  potential  to  alter  the  reactivity  of  the  iron  cell.  Iron  surfaces 
may  become  less  accessible  to  the  aqueous  components  in  the  groundwater,  and/or 
the  precipitated  secondary  minerals  may  provide  new  sites  for  sorption  and  nucle- 
ation.  The  modeling  results  imderscore  the  need  for  more  detailed  studies  of  the 
nucleation  mechanics  of  secondary  minerals  on  Fe“  surfaces. 

Ammonia  was  not  part  of  the  sampling  program  but  was  included  in  the 
modeling  because  it  is  a  hazardous  chemical  and  part  of  the  nitrate  redox  couple 
that  was  postulated  to  be  active  in  the  iron  cell.  The  modeling  shows  ammonia  to 
be  a  thermodynamically  favored  component  that  is  consistent  with  the  disappear¬ 
ance  of  nitrate  in  the  iron  cell  and  the  lack  of  significant  change  to  the  N2(aq) 
concentrations.  It  should  be  considered  for  field  sampling  analyses  at  iron  barrier 
sites  whenever  nitrate  is  known  to  be  present. 

Model  assessment.  The  reactive  transport  system  is  the  amalgam  of  several 
interdependent  processes.  The  ability  to  mechanistically  simulate  the  multifaceted 
behavior  of  this  system  without  significant  calibration  reflects  the  quality  and 
consistency  of  the  data,  and  builds  confidence  in  the  conceptualization  of  the 
system. 

Modeling  provided  a  way  to  test  the  understanding  of  reactive  transport  at  the 
Moffett  Field  Site  and  a  way  to  systematically  identify  plausible  process  mecha¬ 
nisms.  The  flow  through  the  reactive  barrier  and  the  selection  of  chemical  compo¬ 
nents  are  based  on  field  observations.  Once  these  decisions  are  made,  there  are 
very  few  “knobs”  to  tune  the  model.  The  stability  constants  for  all  the  reactions 
are  from  well-established  databases  and  are  not  adjustable.  Furthermore,  boxmdary 
and  initial  concentrations  for  the  selected  species  are  constrained  by  total  concen¬ 
trations  measured  in  the  field  and  speciation  that  satisfy  the  reaction  thermo¬ 
dynamics.  In  this  study,  the  reaction  rates  are  derived  firom  the  literature,  not 
calibrated.  The  fidelity  of  the  model  predictions  with  field  observations  is  good; 
and  where  discrepancies  do  exist,  the  mechanistic  modeling  approach  provides 
insight  on  what  processes  need  to  be  examined  and  how  they  might  be  updated. 

Using  36  chemical  components  and  8  minerals  to  model  a  system  may  appear 
to  be  too  complicated  for  routine  engineering  of  an  in  situ  reactive  barrier;  how¬ 
ever,  the  model  was  built  up  in  a  straightforward  fashion  from  the  observed 
behavior  in  the  system.  Because  of  the  systematic  coupling  of  transport  and  reac¬ 
tions  for  multiple  chemical  components,  simply  adjusting  precipitation  rates  to 
match  the  disappearance  of  certain  ions  will  not  necessarily  be  consistent  with 
other  measured  properties  of  the  reactive  cell  (e.g.,  pH,  alkalinity).  Furthermore, 
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the  reaction  rates  are  operating  within  the  context  of  flowing  groundwater.  The 
character  of  spatial  variations  in  species  concentrations  in  the  reactive  cell  is  the 
direct  result  of  the  interplay  between  rates  of  transport  and  rates  of  reaction. 

Engineered  in  situ  reactive  cells  are  highly  amenable  to  reactive  transport 
modeling.  The  system  being  modeled  is  better  controlled  and  characterized  than 
most  environmental  situations.  In  this  case,  the  chemical  components  and  their 
concentrations  entering  the  reactive  cell  were  comprehensively  studied  and  Fe^ 
was  initially  the  only  mineral  in  the  cell. 


GMS  Implementation 

In  the  GMS,  the  OS3D  reactive  transport  simulator  can  be  used  in  a  stand¬ 
alone  mode  or  in  conjunction  with  multidimensional  flow  fields  generated  by  the 
MODFLOW  groundwater  flow  simulator.  In  standalone  mode,  OS3D  can  solve 
batch  and  reaction  path  chemistry  and  reactive  transport  systems  with  steady, 
uniform  velocity  fields  in  one,  two,  and  three  dimensions.  OS3D  coupling  with 
the  MODFLOW  groundwater  flow  simulation  application  allows  for  reactive 
transport  in  steady  or  transient,  nonuniform,  multidimensional  flow  fields.  This 
linkage  occurs  through  (a)  the  transfer  of  hydrologic  problem  specification 
parameters  from  MODFLOW  (e.g.,  model  dimensionality  and  geometry,  time¬ 
stepping,  spatial  distribution  of  properties)  to  OS3D  using  the  GMS  interface;  and 
(b)  the  generation  of  output  files  (velocities,  sources,  sinks)  by  the  GMS  version 
of  MODFLOW  that  are  read  by  OS3D. 

The  GMS  interface  for  OS3D  was  still  under  development  at  the  time  this 
report  was  written,  but  the  plan  is  to  have  the  GMS  interface  take  the  user  through 
the  OS3D  input  data  specification,  which  includes  solution  options  as  well  as  the 
transport  and  chemistry.  Default  values  will  exist  for  all  solution  options.  Gen¬ 
erally,  the  principal  tasks  for  the  user  will  be  to  specify  the  transport  parameters, 
primary  and  secondary  chemical  species,  minerals,  kinetic  rates,  and  constraints 
for  the  geochemical  conditions.  Transport  parameters  that  can  be  specified  by  the 
user  will  be  spatially  variable  diffusion  coefficients,  dispersivities,  and  the  ability 
to  account  for  tortuosity  through  formation  factors  as  well  as  cementation  factors. 
These  can  be  specified  using  the  GMS  interface  for  cell  properties.  Temperature- 
dependent  diffusion  can  be  specified  with  an  activation  energy. 

For  the  chemistry,  the  user  will  be  presented  with  a  list  of  primary  species  to 
select  from.  Once  the  list  of  primary  species  is  set,  the  user  will  then  be  presented 
with  a  list  of  all  possible  secondary  species  to  select  fi-om.  The  user  has  the  option 
to  select  all  the  possible  secondary  species,  which  can  be  useful  during  the  pre¬ 
liminary  screening  of  species  for  inclusion  in  the  model.  With  the  primary  and 
secondary  species  lists,  OS3D  will  automatically  define  the  stoichiometry  and 
stability  constants  for  each  reaction  (including  minerals).  The  user  can  then  build 
speciated  geochemical  conditions  (i.e.,  sets  of  equilibrated  species  concentrations) 
that  honor  particular  constraints  for  the  primary  species.  For  each  geochemical 
condition,  the  GMS  interface  will  present  the  user  with  a  choice  of  one  of  sk 
possible  constraints  that  are  to  be  used  for  each  of  the  primary  species:  fixed  total 
concentration,  fixed  species  concentration,  fixed  activity,  equilibrium  with  a 
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selected  mineral,  equilibrium  with  a  selected  gas  at  a  selected  partial  pressure,  and 
concentration  that  charge-balances  the  solution.  It  is  at  this  point  that  the  user  will 
eliminate  species  with  relatively  small  concentrations  from  the  final  set  of 
modeled  species.  The  GMS  interface  will  allow  the  user  to  iteratively  build, 
modify,  and  speciate  a  set  of  chemical  components  for  the  modeling  study.  Once 
the  user  is  satisfied  with  the  set  of  chemical  components  and  speciated  geochem¬ 
ical  conditions,  the  conditions  can  then  be  assigned  to  specific  cell  locations  as 
boundary  and/or  initial  conditions. 

A  restart  option  can  be  invoked  that  will  allow  a  previously  run  simulation  to 
continue  from  the  last  saved  time  plane  when  it  was  originally  terminated.  The 
saved  restart  file  can  also  be  used  as  the  initial  condition  for  a  new  sunulation  that 
uses  the  same  problem  geometry  and  chemical  components.  Results  that  can  be 
reported  include  total  and  individual  species  concentrations;  mineral  saturation 
state,  reaction  rates,  and  volume  fractions;  and  aqueous  reaction  rates.  A  variety  of 
options  for  writing  results  to  files  are  available  to  the  user  including  time  series  of 
selected  component  concentrations  at  specific  locations  and  snapshots  of  the 
entire  spatial  domain  at  a  given  time.  GMS  visualization  utilities  can  then  be  used 
to  examine  spatial  and  temporal  distributions  of  selected  parameters. 


Recommendations 

OS3D  is  recommended  for  use  on  multidimensional  advection-diffusion- 
reaction  problems  where  the  numerical  stability  criteria  for  the  advective  trans¬ 
port  scheme  does  not  severely  limit  the  time-step  size  relative  to  the  desired 
simulation  periods.  It  appears  to  be  an  ideal  design  and  performance  assessment 
tool  for  in  situ  chemical  treatment  as  well  as  other  subsurface  characterization  and 
engineering  activities. 

The  OS3D  simulator  can  be  operated  by  users  of  varying  expertise  in  reactive 
transport  modeling  but  does  require  familiarity  with  basic  concepts  of  hydro¬ 
geology  and  aqueous  chemistry.  The  user  is  expected  to  specify  reasonable  and 
appropriate  sets  of  primary  and  secondary  species,  minerals,  reaction  kinetics,  and 
sorption  parameters,  as  well  as  meaningfiil  geochemical  conditions. 

The  user  should  understand  that  while  the  thermodynamics  of  most  inorganic 
chemical  reactions  of  interest  are  well  known  and  exist  in  the  EQ3/EQ6  database, 
some  of  these  reactions  are  kinetically  controlled  with  rate  laws  that  are  generally 
developed  xmder  specific  experimental  conditions.  To  get  the  most  benefit  from 
the  OS3D  simulation  technology,  the  user  should  become  familiar  with  the  sup¬ 
ported  rate  law  forms  and  be  prepared  to  modify  and  augment  the  default 
database. 

The  solver  technology  in  OS3D  has  been  shown  to  be  quite  robust  in  solving 
a  variety  of  nonlinear  reaction  systems  under  computationally  demanding  condi¬ 
tions;  however,  convergence  is  not  guaranteed.  It  is  recommend  that  users  not 
adjust  the  convergence  criteria  unless  they  have  a  full  appreciation  for  the  poten¬ 
tial  consequences  these  adjustments  may  have  on  the  accuracy  of  the  simulation. 


Chapter  5  OS3D 


65 


Because  of  the  large  number  of  chemical  components  that  might  be  employed, 
a  minimum  of  128  MB  of  memory  for  general  use  and  256-512  MB  of  memory 
for  larger  three-dimensional  problems  are  recommended.  Large,  high-resolution 
three-dimensional  problems  should  be  run  using  the  parallel  processing  version  of 


OS3D. 
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6  Conclusions 


Four  simulators  for  modeliug  in  situ  subsmface  remediation  were  advanced 
and  incorporated  into  the  GMS  during  this  project:  MT3DMS/SEAM3D, 
UTCHEM,  NUFT,  and  OS3D.  The  latest  version  of  the  OS3D  code  is  now  being 
referred  to  as  CRUNCH.  These  four  simulators  are  a  significant  advancement  in 
the  present  status  of  subsurface  modeling  for  remediation  evaluation  and  design. 
Each  of  these  four  simulators  has  unique  capabilities  for  simulating  innovative  in 
situ  treatment  technologies.  There  is  some  overlap  of  these  models  in  terms  of 
capabilities,  but  this  is  viewed  as  an  advantage. 

The  MT3DMS  code  is  a  multispecies  transport  code  with  significant  advance¬ 
ments  over  the  previous  version,  MT3D.  MT3DMS  is  the  transport  module  for 
SEAM3D  and  RT3D,  both  reactive  models  in  the  GMS.  SEAM3D  has  many 
features  for  simulating  intrinsic  and  engineered  bioremediation.  The  model  can 
handle  the  full  suite  of  aerobic  and  anaerobic  degradation  processes  for  petroleiim 
and  chlorinated  hydrocarbons  under  saturated  flow  conditions. 

The  UTCHEM  code  is  probably  the  most  versatile  of  the  simulators  as  it 
allows  for  multiphase,  multicomponent  flow  and  transport.  The  primary  emphasis 
of  this  model  is  simulating  surfactant  and  cosolvent  flushing  of  the  oil  phase. 
UTCHEM  is  not  well  suited  for  simulating  some  other  multiphase  flow  problems, 
such  as  soil  venting  and  steam  injection.  However,  the  model  also  has  modules  for 
simulating  bioremediation  and  geochemical  reactions,  although  features  for  the 
latter  module  have  not  been  interfaced  with  the  GMS.  UTCHEM  will  handle 
variably  saturated  flow  conditions. 

The  NUFT  code  also  simulates  multiphase,  multicomponent  flow  and  trans¬ 
port.  NUFT  was  originally  developed  for  simulating  steam  injection  and  soil 
heating  and  venting  technologies  for  in  situ  remediation.  Thus,  this  code  is  well 
suited  for  these  types  of  applications,  but  it  is  not  well  suited  for  modeling  sur¬ 
factant  and  cosolvent  flushing  technologies  for  oil  removal.  NUFT  also  has  a  bio¬ 
reactions  component,  but  this  feature  has  not  been  interfaced  -with  the  GMS. 
NUFT  will  handle  variably  saturated  flow  conditions. 

The  OS3D  code  is  a  versatile  geochemical  code  best  suited  for  in  situ  chemi¬ 
cal  treatment  technologies,  such  as  permeable  treatment  walls.  This  model  is 
intended  only  for  application  to  groundwater  plumes  under  saturated  flow 
conditions. 
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At  the  time  of  this  publication,  three  of  the  four  featured  simulators,  plus 
MT3DMS,  have  been  implemented  into  the  GMS.  Work  on  implementing  OS3D 
was  still  in  progress.  Tutorials  for  all  of  these  simulators  either  have  been  or  are 
being  developed  to  aid  users  in  applications.  As  a  result  of  this  project,  the  GMS 
now  has  extensive  capabilities  for  modeling  a  wide  array  of  subsurface  and 
grormdwater  flow  and  reactive  transport  problems  involving  a  host  of  cleanup 
and/or  management  alternatives. 

Some  of  the  advanced  features  of  these  simulators  will  present  a  challenge  for 
less  experienced  groimdwater  modelers.  Thus,  it  is  recommended  that  modelers 
take  the  time  to  become  closely  acquainted  with  the  model  processes  being  used. 
Due  to  the  comprehensive  features  of  these  four  simulators,  these  models  are  well 
suited  for  detailed  design  and  evaluation  of  in  situ  treatment  alternatives  prior  to 
field-scale  implementation.  However,  with  the  Map  Tool  and  other  site 
conceptualization/characterization  features  of  the  GMS,  once  a  site  study  has  been 
set  up,  it  is  relatively  straightforward  to  apply  various  GMS  models  to  the  site, 
enabling  multiple  simulation  alternatives  to  be  conducted  and  evaluated  expedi¬ 
tiously,  even  diuing  the  alternatives  screening  phase  of  a  project.  The  GMS 
models  will  provide  an  invaluable  tool  for  improving  subsurface,  in  situ  remedi¬ 
ation  selection,  design,  and  operation,  thus  potentially  providing  significant 
reductions  in  site  cleanup  costs. 
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